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Abstract 

We review several facets of the hydrodynamic description of the relativistic heavy ion collisions, 
starting from the historical motivation to the present understandings of the observed collective 
aspects of experimental data, especially those of the most recent RHIC and LHC results. In 
this report, we particularly focus on the conceptual questions and the physical foundations of the 
validity of the hydrodynamic approach itself. We also discuss recent efforts to clarify some of 
the points in this direction, such as the various forms of derivations of relativistic hydrodynamics 
together with the limitations intrinsic to the traditional approaches, variational approaches, known 
analytic solutions for special cases, and several new theoretical developments. Throughout this 
review, we stress the role of course-graining procedure in the hydrodynamic description and discuss 
its relation to the physical observables through the analysis of a hydrodynamic mapping of a 
microscopic transport model. Several questions to be answered to clarify the physics of collective 
phenomena in the relativistic heavy ion collisions are pointed out. 


1 Introduction 

Hydrodynamics is a theoretical framework to describe the motion of fluids based on their local property 
and the conservation laws of energy and momentum, and other conserved quantities, for example mass 
in non-relativistic hydrodynamics. The main advantage of hydrodynamics resides in the fact that a huge 
number of degrees of freedom contained in the microscopic composition of the fluids is drastically reduced 
to a few macroscopic hydrodynamic variables which represent the local property of the fluid. We refer 
to this property of hydrodynamic description as “locality”. The local property is usually represented 
by thermodynamic relations among the hydrodynamic variables frequently called as equation of state 
(EoS), plus transport coefficients. That is, it is assumed that the fluid is in the state of local thermal 
equilibrium (LTE), or very close to it (see Sec. [2]). 

Hydrodynamics is one of the oldest classical phenomenological theories which has played an impor¬ 
tant role in the development of science and technology. Its relativistic form was already developed in 
the early stages of the theory of relativity DE] mainly in the context of astrophysics and cosmology. As 
we will report in this review, applications of hydrodynamics to the microscopic system such as nuclear 
collective motion and heavy ion collisions have been shown to be successful. It is amazing to observe 
that such a simple scheme of dynamics is able to cover various phenomena from cosmological to hadronic 
scales. 

When we think of typical dynamics of a fluid, we immediately imagine its flow profile as we usually 
observe in daily life, such like waves of sea, ripples of wine in a glass, vortices in a water-sink, turbulences 
in winds, etc. Indeed, in most of these phenomena, the hydrodynamic picture is known to be well 
established since the validity of the concept of locality in describing the properties of fluids through 
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LTE seems obvious due to the huge number (~10 20 ) of particles involved in a visible dynamics of these 
fluids. In such cases, if the property of the matter is given we can study its response as an initial 
value problem. Sometimes we can use the hydrodynamic description to determine the initial condition 
which leads to the specified final state after the hydrodynamic evolution. On the other hand, if the 
properties of the matter is not known, then we can introduce a model to represent them in terms 
of a few hydrodynamic parameters, and infer them by comparing the predictions with experiments, 
provided that the hydrodynamic picture is in fact valid. Of course, hydrodynamic responses of the 
matter in principle can be nonlinear, and sometimes they are extremely complex. In such cases the 
above mentioned hydrodynamic modeling to deduce the initial condition or properties of the matter in 
question is not trivial at all. In particular, for example in the presence of turbulence, the time evolution 
becomes chaotic and the one-to-one correspondence of the initial and final states is lost [3]. 

As mentioned, the hydrodynamic picture is applicable for a vast class of different scales. However, 
it should be kept in mind that the meaning of the locality changes depending on the scale of the 
system in question. For example, if we try to study hydrodynamics of the atmosphere for weather 
forecasting purpose, larger-scale simulations are implemented to cover the significant area of the earth. 
There, the locality means that any inhomogeneity of the air is completely neglected within the volume, 
say, of the order of m 3 , and in fact, too precise information is not required. However, this does not 
necessarily mean that small scale dynamics is completely neglected. For example, in the presence of 
turbulences in a smaller volume than the required precision, they are smeared out and counted as the 
internal degrees of freedom of the air in our observational scale. In these situations, the transport 
coefficients might be modified and some effective quantities, e.g. the so-called eddy viscosity should be 
used [3j. Such a situation occurs frequently when we perform numerical simulations where a certain 
scale should eventually be introduced due to discretization of space as mentioned above in the example 
of large scale numerical simulations with cutoff scales, which reproduce the observable phenomena. 
This happens also in astrophysical applications such as supernova, galaxy formation and cosmological 
problems. This means that the hydrodynamic parameters determined from microscopic theories or 
laboratory experiments are not necessarily the same as those used in real simulations. 

It is more than decades ago since the hydrodynamic modeling in relativistic heavy ion collisions has 
been shown to be very successful in understanding the nature of so-called collective flow parameters. 
In particular, after establishing the almost ideal fluid scenario and the picture of the strongly coupled 
deconfinement phase of quarks and gluons (quark and gluon plasma - QGP), the hydrodynamic de¬ 
scription is considered to be indispensable for some stages of the dynamics of the expanding produced 
matter in the relativistic heavy ion collisions. These successes lead to the expectation that some im¬ 
portant bulk properties of the QGP and its initial state just after the collisions can be determined to a 
quantitative level using the hydrodynamic analysis. 

However, the situation of relativistic heavy ion collisions is quite different from the aforementioned 
macroscopic phenomena. First, we are unable to observe directly the time evolution of flow profile of 
the matter produced by the collisions in event-by-event basis (EbyE) but we only observe free-streaming 
final state particles, instead. In fact, the two scenarios, the hydrodynamic flow and the collection of 
free-streaming particles are essentially different. The transition from one scenario to the other is a very 
complex transport process and referred to as “freeze-out”. Second, in the study of relativistic heavy ion 
physics, the initial condition and the property of the matter are exactly what we want to determine. 
In the hydrodynamic approach, these are inputs and are not observable directly determined from the 
experimental data. Thus we have to resort to theoretical models to specify the initial condition and 
the property of the matter, which are out of the scope of the hydrodynamic formalism. Finally, at 
the present stages, we are not able to be sure on why and when a hydrodynamic description becomes 
valid in such an extremely explosive dynamics in the scale of subatomic systems. All of the above three 
aspects of the hydrodynamic approach constitute very challenging problems, and extensive studies have 
been done both experimentally and theoretically. 
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For the sake of bookkeeping, in this report, we classify these studies roughly in the following three 
categories. 

1) Assume the validity of hydrodynamic scenario as working hypothesis, construct a unified descrip¬ 
tion of relativistic heavy ion collision, with the help of certain models of initial conditions and 
also the freeze-out processes. By applying such a description to the analysis of experimental data, 
try to determine quantitatively the parameters in the model. 

2) Focus mainly in qualitative aspects of the hydrodynamic picture. Look for robust and less model- 
dependent hydrodynamic signals to explore the initial state dynamics of the collision. Also, 
complementary to the approaches classified in 1), study analytic or simplified solutions to clarify 
the role of hydrodynamic signals. 

3) Establish the foundation and validity of the hydrodynamic approach in the scenario of relativistic 
heavy ion collisions, and also look for possible alternative (or variant) interpretation of collective 
phenomena. 

Of course the above three categories are complementary and also interrelated so that they are 
often not quite separable. The category 1) is presently the main stream of hydrodynamic studies of 
relativistic heavy ion collisions. In fact, it is the most urgent and important task to push forward the 
known hydrodynamic approaches and establish the method of analysis of newly coming experimental 
data. 

On the other hand, this approach deals with the precise quantitative comparison to experimental 
data so that it necessarily involves quite complex modelings for the initial condition and for the final 
state interactions, in addition to rather sophisticated numerical techniques. Thus sometimes the effects 
of physical parameters in the model are not quite evident for those who do not have these resources 
at their disposal to study the individual technical issues. Furthermore, the space of parameters to be 
determined is quite large so that the uniqueness of the solution is not guaranteed. In this aspect, the 
approaches of the category 2) and 3) are helpful to characterize the flow dynamics from more general 
insights giving feedback to the approach 1). They are also useful to extend the hydrodynamic picture 
in different situations such as the beam energy-scan program (BES) at RHIC/BNL, FAIR/GSI and 
NIC A/ JINR [U 15] . In addition, through these approaches we can clarify the limitations and applicability 
of relativistic hydrodynamics. This is a fundamental point for the sound understanding of physics of 
relativistic heavy ion collisions, because the reproduction of experimental data is of course not the 
ultimate objective of the hydrodynamic approach. Extensive studies on the role of the hydrodynamic 
approach in relativistic heavy ion collisions have been done, such as these comprehensive reviews and 
books [SHU]. In particular, these recently published reviews HEHU] report the present status, giving 
an overview of the state-of-art of the hydrodynamic approach mainly focusing on 1), including LHC 
results. Therefore, in this work we will not repeat in detail the works already reported and refer 
readers to these review articles and books as well as the original references therein. Instead, we focus 
complementary studies on hydrodynamics in relativistic heavy ion collisions, mainly those to be classified 
in the categories 2) and 3), with the exception of some key point studies which characterize category 1). 
As mentioned, they are equally important for a deep understanding of the physics of collective signals 
in relativistic heavy ion collisions. 

The present paper is organized as follows. In Sec. [2j we give a brief overview on the historical 
motivation and developments of the hydrodynamic approach in relativistic heavy ion collisions. In Sec. 
[3j we show the basic structure of relativistic hydrodynamics and discuss the physical concepts of the 
approach, both for the ideal and dissipative cases. In Sec.[4j we discuss how the hydrodynamic approach 
is applied in practice to relativistic heavy ion collisions. There, we focus on several specific physical 
problems necessary to connect the hydrodynamic variables to the physical observables of relativistic 
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heavy ion collisions, such as the initial condition, equation of state and transport coefficients and the 
freeze-out process. In Sec. [5] we give an overview on how these questions are dealt with in the present 
studies and give a concise summary on the present status and accomplishments of the hydrodynamic 
approach in heavy ion collisions. In Sec. [6j we report on the theoretical efforts from a more general point 
of view, for example, to formulate the dissipative relativistic hydrodynamics from the various methods, 
to find the analytic classes of solutions, to establish foundations of the approach in the scenario of 
relativistic heavy ion collisions, inclusion of some new dynamical effects, etc. In Sec. [7] we focus our 
attention on the aspect of coarse-graining in physical observables and investigate its effects on the origin 
of hydrodynamic description taking an example of event generator code based on a transport theory. 
Sec. [8] is dedicated on the discussion of future perspectives of the hydrodynamic approach in relativistic 
heavy ion collisions. 


2 Relativistic Heavy Ion Collisions and Hydrodynamics 

In this section, we briefly review the historical background of the studies of relativistic heavy ion 
collisions. Research in this area started with theory-driven motivations to search for new, theoretically 
predicted states of matter. Since then, it has developed into an experimentally driven science that works 
not only just to confirm these new states of matter but also to quantitatively determine the properties 
of these states, even leading to new thought provoking questions. 

2.1 Early Days 

The use of thermodynamical concepts in high energy elementary particle collisions was introduced for 
the first time by Fermi to model the multi-particle productions [21]. Landau extended Fermi’s fireball 
model to analyze the effects of collective motion of the produced particles in the framework of relativistic 
hydrodynamics [221 23] • He showed the importance of the longitudinal expansion of the fireball using 
the Khalatnikov’s analytic solution of (1+1)D relativistic hydrodynamics [25]. At that time, of course, 
there was no idea of the mechanism of high energy hadronic collisions which leads to multiple pion 
productions. Thus, instead of modeling the production mechanism, Landau simply assumed that the 
created matter could be regarded as a hot gas of massless pions in thermal equilibrium as Fermi did 
(which is not much far from the present idea of the parton gas). But he argued that the fireball created 
in a proton-proton (p-p) collision suffers a strong longitudinal, rather than isotropic, expansion due to 
the large pressure gradient in this direction because of the Lorentz contraction in the initial state. This 
affects naturally the energy dependence of the multiplicity of final particles on the incident energy and 
generates the characteristic rapidity distribution. 

Fermi and Landau’s thermodynamic/hydrodynamic approaches have been pursued by many authors 
to understand the nature of multiple particle production in hadronic collisions observed in cosmic ray 
phenomena. These works served as the conceptual origin in the studies of bulk aspects of the strongly 
interacting matter at very high energies. In the middle of 60’s Hagedorn predicted the existence of the 
limiting temperature for the hadron resonance gas (HRG) [25], and later this behavior was discussed 
with the deconfinement of quarks |26j (see also Ref. G3). Recently, his statistical bootstrap approach 
was revisited from the modern perspective of QGP to investigate the hadronization mechanism and 
transport properties of HRG [I2HMSU] . Also a theoretical framework of HRG was developed in Ref. pTi] • 
In his hydrodynamic studies on multi-particle production in high energy hadronic collisions, Carruthers 
pointed out in 1974 [32] that the matter which obeys hydrodynamics does not need to be a fluid of known 
particles. Instead, he emphasized the importance of using the coarse-grained macroscopic quantities 
such as energy-momentum tensor and conserved densities to describe the possible new states of matter, 
which he calls pre-matter. At that time, the concept of freeze-out to derive the observable particles 
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spectra from the hydrodynamic picture was already studied and the well-known Cooper-Frye formula 
was developed in this context for the multiple particle productions in high energy hadronic collisions 


(see the discussion in Sec. 4.3). 

The hydrodynamic approach has also been used in nuclear physics since the early days, although 
before the 70’s the main interests have been concentrated on non-relativistic energies. This is because 
the principal object of nuclear physics at that time was to investigate the properties of nuclei and their 
reactions from the point of view of many-body quantum systems. The underlying nuclear Hamiltonian 
is given by two-body nucleon-nucleon interactions, represented by potentials which were not established. 
Therefore, for the nuclear physics community, scientific interests on highly compressed or excited nuclear 
matter were not popular yet. On the other hand, for the elementary particle physics community, high 
energy proton-nucleus (p-A) or nucleus-nucleus (A-A) collisions were considered to be too complicated, 
since they certainly provoke excitations of many degrees of freedom and the study of yet unknown 
hadronic interactions through such complex processes seemed to be almost hopeless. 

Nevertheless several pioneering studies emerged suggesting possible existence of exotic states in 
highly excited hadronic gas as mentioned above, or very high density nuclear isomer states [34]. The 
use of high energy A-A collisions to investigate the nuclear EoS was also discussed [35, [36] . In the early 
70’s, the first relativistic nuclear beam was realized at LBL 03- 

At the same time, in 1968, the discovery of the Erst pulsar [38, [39] and its identification as a 
neutron star, which was predicted long before by Landau na as the form of giant nucleus bound by its 
gravitational force (more precisely postulated by Baade and Zwicky |4I] after the discovery of neutrons), 
provoked strong interest in the properties of highly compressed nuclear matter in the nuclear physics 
community (see for example Ref. [12] )■ The possibility of a quark matter star was first suggested 
in 1970 by N. Itoh [43]. The program of relativistic heavy ion collisions to study the properties of 
matter at extreme conditions was advocated by T.D. Lee in 1974, which is now grown up to one of the 
contemporary frontiers of sciences. For more details of these developments, we refer readers to Refs. [44] 
and [45]. It should be noticed that all these happened at the time when Quantum Chromodynamics 
(QCD) was emerging [46, 47] (see also Ref. [08] for the development of QCD) and asymptotic freedom 
was just discovered m [50] . but well before the experimental discovery of jets [51] [52]. Of course 
now the well-known concept of Quark and Gluon Plasma (QGP) [53j was even not existing among the 
community. 

Thus, in the 70’s and 80’s, several particle accelerators were adopted to accelerate nuclear projectiles, 
e.g, BEVALAC/LBL, DUBNA, AGS/BNL, SPS/CERN. Initially the incident energy of the projectiles 
for fixed target was around 1 GeV per nucleon, and to investigate whether such a state of the compressed 
nuclear matter can be achieved was one of the main topics. During this time period, several relativistic 
effects and methods were developed, such as relativistic mean held theory and relativistic intranuclear 
cascades [Mli55], in addition to the application of relativistic hydrodynamics. It was first pointed out by 
Stocker et al. [56], that the side peaked angular distribution of protons favors the hydrodynamic scenario 
compared to a simple binary intranuclear cascade. With the theoretical developments of QCD and the 
discovery of asymptotic freedom, partons were identified with quarks and gluons in this regime. Thus, 
the concept of a plasma of quarks and gluons (QGP) became more concrete as a new state of matter 
at extreme conditions. Many theoretical investigations on possible QGP signals in relativistic heavy 
ion collisions were proposed, such as strangeness enhancement [57], J/T suppression [58], restoration 
of chiral symmetry [52], jet quenching j6UJ and so on mm- Influences of the possible new state of 
matter on collective how dynamics were also discussed [62, 64]. For more details and further reading of 
the early stage of the relativistic heavy ion collisions, see Refs. [651167] and books [68] [69] which convey 
the excitements and perspectives for the new research area in nuclear physics of these times. 

In order to have a clear signal predicted from theories, we need to achieve a state of QGP in thermal 
equilibrium, or in other words, we need a large space-time scale of reaction processes. This is why we 
consider heavy ion collisions the most promising area to study QGP in the lab. Elementary collisions 
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Figure 1: (color online) QCD Phase Diagram conjectured from QCD. The left and right panels represents 
the two distinct predictions, with and without the presence of critical points, respectively. Figur^Jtaken 
from Ref. [72] • 


such as p-p are considered to have space-time scales too small to reach thermal equilibrium in created 
matter (for recent surprises, see Sec. [5]). ffeavy ion collisions play the role of a “pressure cooker” to 
maintain matter in a hot and dense enough state to produce QGP. In contrast to elementary collisions, 
the geometry of the initial configuration plays an essential role in heavy-ion collisions. Of course, the 
initial collision geometry is not a direct observable but can be inferred from the final state multiplicity of 
produced particles. Nowadays, developments of experimental techniques and findings of new methods 
and observables make it possible to determine the centrality class of the collision in terms of global event 
observables, such as particle multiplicity and forward energy measured by the zero degree calorimeter 
(ZDC). Also techniques of particle identification and methods of particle correlation measurements 
allow us to obtain more exclusive data, which enables precise comparisons with theoretical models at 
the quantitative level. The recent review papers [701 Ej by the authors working directly on collective 
flow phenomena in relativistic heavy ion collisions summarize these developments through the last four 
decades, from the very early days to the present. 

2.2 QCD Phase Diagram 

As mentioned above, one of the main objectives of the relativistic heavy ion program is to investigate 
the properties and dynamics of the QCD matter at extreme conditions of temperature and density. The 
present image of the phase diagram of QCD is summarized in Fig. [l] 

As seen in Fig. [lj the QCD phase diagram has a much more complex structure compared to the 
early speculations deduced from the bag model. First, the lattice QCD (1QCD) calculations showed that 
the transition from the hadronic phase to the QGP phase at the vanishing baryonic chemical potential 
(Hb = 0) is a smooth cross-over and not a first order transition [73]. This behavior is now considered to 
be well established and we will come back to this point later with respect to the QCD EoS. Second, the 
order of the phase transition depends on the baryon chemical potential, leading to a possible critical 
point (left panel in Fig. [Tj) . The existence of this kind of critical point was first suggested in 1989 
m in the context of chiral symmetry breaking. Although some finite chemical potential results also 
show the existence of a critical point [75, 76], this is still an open ended question m, as shown in the 
right panel of Fig. [I] There, the critical point does not even exist anywhere m- Third, the chiral 
phase transition does not necessarily go with the confinement phase transition. Possibilities of different 
behavior of the two phase transitions have been discussed in 1QCD and in some theoretical models, and 
the possible existence of a new phase called quarkyonic matter was suggested in Ref. [75] where the 

“Reprinted from Prog. Part. Nucl. Phys. 72, K. Fukushima and C. Sasaki, The phase diagram of nuclear and quark 
matter at high baryon density, 99, Copyright © 2013, with permission from Elsevier. 
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Figure 2: Thermal Model Fit of data from Pb-Pb at a/snn = 2-76 TeV. Figure^] taken from Ref. [90]. 

confined state with the chiral symmetry is realized. For more detailed information, see review papers 
HUES]. For lower temperatures and high baryonic chemical potentials, quark matter as conjectured 
by Itoh [43] appears, but still exhibits complicated structures such as various color-superconducting 
phases. For more detailed discussions on the theoretical investigations of the QCD phase diagram, see 
the recent reviews, Refs. [721 ng. See also Ref. [RT] for a general review on statistical properties of 
strongly interacting matter. 

One of the basic motivations of relativistic heavy ion collisions is to experimentally investigate 
the presumed rich structure of the QCD phase diagram. By definition, the phase diagram represents 
the states of a homogeneous matter in thermal equilibrium with infinite volume. As mentioned, it is 
not obvious whether the produced matter achieves a thermally equilibrated state in relativistic heavy 
ion collisions. One of the experimental indications of thermal equilibrium is observed in the relative 
abundances of produced particles. In the thermal models, the relative abundances of various hadron 
species are given by the thermal equilibrium distribution at the so-called chemical freeze-out stage, 
specified only by the two parameters, the temperature T and the baryonic chemical potential hb [82ti88] . 
In more sophisticated models, effects of the chemical non-equilibrium are also considered by introducing 
fugacities for quarks [52]- In Fig. [2] we illustrate the recent results of particle ratios with the typical 
thermal model fits for the RHIC and LHC energies. 

Whether the perfect chemical equilibrium is achieved or not, it is a fact that quite good reproductions 
of the data of particle ratios are obtained in terms of the thermal approach. This indicates that the 
chemical freeze-out process occurs as if the system is in a state very close to the thermal equilibrium in 
relativistic heavy ion collisions. On the other hand, it is intriguing that the thermal models are shown 
to work also well for the particle ratios in high energy p-p and e-e reactions [E1E2]- Note that in these 
examples, the thermodynamic model is applied to the system averaged over many events. 

The thermal model above is applied to the final HRG. Thus, this does not imply directly that the 
system reaches thermal equilibrium at the QGP stage. Furthermore, if the particle ratios reflect the 
complete statistical equilibrium of a pure hadron gas, then the information of QGP should be washed 
out. The signals of QGP in flavor enhancement/suppression such as strangeness enhancement or J/T 
suppression should stand out from this thermal model fit of the pure HRG. See more detailed discussion, 
Chap. 24 of Ref. [13] and references therein. 

“Reprinted from J. Stachel et al., Confronting LHC data with the statistical hadronization model , J. Phys. Conf. 
Ser. 509 (2014) 012019, doi:10.1088/1742-6596/509/1/012019, under the terms of the Creative Commons Attribution 3.0 
license. 
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One of the QGP signals which may survive throughout final state hadron interactions is collective 
flow. Loosely speaking, if the QGP phase emerges, the corresponding pressure is much larger than 
that of a hadron gas for a given energy density, so that in the presence of QGP phase any spatial 
inhomogeneity enhances an acceleration of the expansion of the matter. Such an acceleration should 
manifest in the final particle angular distribution. 

In order to quantify the above expectation, we have to model the collision process in terms of a 
hydrodynamic description. In the next section, we show the basic equations of relativistic hydrodynam¬ 
ics in the most concise form as possible. More detailed physical aspects involved in the hydrodynamic 
formulation will be discussed in the following sections. 

3 Structure of Relativistic Hydrodynamics 

In this section we introduce the basic concepts and equations of relativistic hydrodynamics. For the 
sake of a pedagogical purpose, we start from the familiar non-relativistic Euler equation for the ideal 
fluid to clarify the concept of coarse-graining used in the hydrodynamic description. 

3.1 Non-relativistic Euler equation 

The non-relativistic Euler equation for an ideal fluid has the following form, 

v + (v-V) v = -—VP, (1) 

Pm 

where v is the velocity field, p rn is the mass density and P is the pressure. In basic physics text books 
this equation is usually derived by applying Newton’s equation of motion to a small cell of the fluid and 
using Pascal’s law which claims that the pressure is isotropic. Here, we derive it from the conservation 
of momentum in order to clarify the physical structure of the relativistic hydrodynamics. To do this, we 
have to introduce several physical assumptions: first, states of the matter, independent of the details 
of internal states of microscopic degrees of freedom, can approximately be described in terms of a 
few number of macroscopic classical fields. We refer to those as hydrodynamical variables, and they 
are constructed at an arbitrary space-time point by taking the averages over the microscopic degrees of 
freedom within the vicinity of the point. The size of vicinity should be specified by a certain scale, called 
the scale of coarse-graining (C-G). We thus represent the system as a continuum medium, but remind 
that for each space-point we have always to associate a small volume determined by C-G size, to which 
we refer as “fluid element”. By definition these fields should be insensitive to a microscopic dynamics 
which occurs inside the fluid element, so that any hydrodynamic variables in the co-moving frame should 
be almost constant, and isotropic in the absence of external forces or internal degrees of freedom such 
as spin. Quantitatively speaking, C-G scale should be less than the minimum of 1/ |Vln(<5)| where Q 
is the typical hydrodynamic variables. The condition of the continuum limit for the case of a gas of 
particles, is achieved when the mean-free path is much less than the C-G scale. The ratio of these two 
numbers is called Knudsen number K. Second, forces exert on one fluid element come only from the 
adjacent fluid elements and expressed in the form of stress tensor. Third, there should be no tangential 
stress for a fluid at rest. This condition is crucial to distinguish the two possible continuum media, 
say, fluid and solid. In particular for the ideal fluid where we discuss here, this condition is valid for 
any fluid element. To constitute a self-contained dynamical system under the above conditions, we 
have to construct a closed system of dynamical equations with these classical fields. This last condition 
implies that the scales for macroscopic dynamics defined by C-G should be clearly separated from the 
microscopic ones. 

The important question is how to choose the hydrodynamic variables. Not every physical quantity, 
even after C-G is necessarily appropriate for the description of dynamics. Normally, conserved densities 



are good candidates. To see this, let n and J be, respectively, some conserved density and its current, 
which are expressed in microscopic degrees of freedom and satisfy the continuity equation. Then the 
corresponding C-G density and current, n and J, should also satisfy the continuity equation, 

d t n = -V • J. (2) 


Suppose the C-G scale is given by L. Then from the first condition we mentioned above, the right hand 
side of Eq. Q is at most the order of |J| /L. In such a case, a significant change of n occurs only for 
the time-scale much larger than that of the microscopic motion. This means that the dynamical time 
scale of n becomes also that of macroscopic scale compatible with the C-G size. Therefore, at least 
n is appropriate variable to characterize the macroscopic motion of the system (see a more detailed 
discussion in Secs. 6.1.2 and 6.1.3). In general, the current J can still contain microscopic motions. 
In the case of the mass conservation in the non-relativistic hydrodynamics, it can be shown that the 
corresponding current is also a macroscopic quantity associated with the momentum conservation. 

In the non-relativistic case, the mass is a conserved quantity and the corresponding continuity 
equation for the mass density p m and the current j = p m v is given by, 


dpm 

~df 


-V ■ (p m v), 


(3) 


where v is the velocity field. For the later purpose, we express this continuity equation in a Lorentz 
covariant form, d^j 11 = 0, where (j M ) = ( p m c p m v ), with ( x = (cf, r) and c is the speed of light. 

The above continuity equation can readily be generalized for conserved vector quantities. For exam¬ 
ple, the total momentum of the medium, f dV p m v should be conserved in the absence of any external 
forces. Thus in terms of the continuity equation, we have 


d t (pmV l ) + d k (p m v l v k ) = -diP, 


(4) 


where the right-hand side accounts for the interaction among the fluid elements through the surfaces 
with P being pressure. With the use of Eq. (|3]), we arrive at Eq. (JTj). 

Equations ([TJ) and (J3| constitute a system of partial differential equations for p m and v but due to 
the presence of the pressure P, the system is not closed yet and we need another condition. For the 
case of a barotropic fluid where P is given by a function of p m , we can solve the system without any 
other knowledge of the matter, such as internal energy. 

To be complete, let us consider the energy conservation of the fluid element. The time derivative of 
the energy can also be expressed as the sum of the energy flux through the surface and the rate of the 
work done by the external forces and other fluid elements. The energy density is given by p m v 2 / 2 + e int , 
where the first term is the translational kinetic energy density and the second term comes from the 
internal energy of the fluid element. The conservation of the energy is then written as 

dt ( pmC 2 + ““V 2 + Eint'j + V ■ V ( p m C 2 + y-V 2 + E int + P j =0, (5) 

where the rest mass energy conservation is added. Finally, we can combine Eq. Q and Eq. ([5]) and 
rewrite them as a “relativistic covariant” expression, 

Wnr = 0. (6) 


where 


rjil^LU _ 

1 NR — 


PmC 2 (1 + /3 2 / 2 + Eint/PmC 2 ) p m C 2 /3 T ^ \ 

PmC 2 (1 + f3 2 /2 + [ £i n t + P] /c 2 ) /3 p m c 2 /3 (3 T + PI J 
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(7) 





and f3 = v/c. Note that the above matrix can be written in the symmetric form as 

_ ( PmC 2 (1 + /3 1 2 /2 + Sint/ PmC 2 ) PmC 2 (1 + /3 2 /2 + \£i n t + P] / p m C 2 ) /3 T 
^ \ PmC 2 (1 + /3 2 /2 + [£j n t + -P] / PmC 2 ) /? PmC 2 f0 (0 T + PI 

without altering the energy and momentum conservation equations in the non-relativistic limit. Equa¬ 
tion which gives the correct Euler equation, looks like a covariant expression. However, differently 
from the conserved charge density in the continuity equation, the above T >w NR is not yet a covariant 
tensor under the Lorentz transformation. The correct covariant quantity is called energy-momentum 
tensor which is defined in the next subsection. 



3.2 Relativistic Euler equation 

Having in mind the physical structure of the non-relativistic case, we derive the relativistic Euler 
equation in a more mathematical way to manifest the Lorentz-covariant structure. We first consider 
the hydrodynamic variables associated with energy-momentum conservation. 

The three-velocity field v now should be extended to the four-velocity u M with the normalization 
condition = 1 . At a given space-time point, we can always find the frame where tA —y ( 1 , 0 , 0 , 0 ). 
We call such a frame local rest frame (LRF). As we saw in the non-relativistic derivation, the energy- 
momentum conservation is described in terms of Eq. So we expect that the relativistic Euler 

equation can be derived by rewriting it in a covariant form. The corresponding quantity, to Tfff 
should be a second rank tensor which transforms under the Lorentz transformation. Furthermore, in 
the ideal case the only non-scalar quantities are the four-velocity field rA and the metric tensor g^, so 
that the tensor structure of this T should be constructed only from these. Then, the most general 
expression of a second rank tensor is written as 


T' w = Au fl u l/ + Bg 


where A and B are scalar functions. In LRF, u M —> (1, 0, 0, 0), we have 


T*" -t- K rf = 


A+B 0 \ 

0 -BI J ' 


(9) 


( 10 ) 


Since we always interpret T 00 as the energy density (including the rest mass and kinetic energy) in any 
reference frame, it should reduce to the proper energy density s in the local rest fram^] 

T° l ° rf = A + B = s. (11) 


Furthermore, comparing Eq. (10) with Eq. 
scalar function 


for 0 — 0 , it is natural to interpret that the proper 


B as the thermodynamic pressure P in LRF. With this identification, we have the 
energy-momentum tensor in a general frame as 




{£ + P) 7 2 - P {£ + P ) 7 2 /3 
{£ + P) 7 2 /3 (e + P ) 7 2 (0/3 t + PI 


{£ + P ) vAu v - Pg' kV . 


( 12 ) 


Since j3 is the (three-)velocity of the fluid element, it determines the Lorentz boost from LRF to the 
observational frame. The associated Lorentz factor is given by 7 = 1/y/l — (3 2 , as usual. The energy 
and momentum conservations are then expressed in terms of this energy-momentum tensor as 


d u T tw = 0. (13) 

1 Proper density of any thermodynamic quantity defined in LFR in relativistic hydrodynamics is by definition an 

invariant quantity. But differently from the non-relativistic case, the corresponding quantity in the observational frame 
has different value. In this paper, to distinguish this observational value from the proper density, we denote it with *. 
For example, the proper charge density in LRF, say n is denoted as n* in the observational frame. 


10 





This equation constitutes four constraints among the five unknown hydrodynamic variables, e, P, and 


j3. If P is expressed as a function of e, then Eq. (13) forms a closed system for four variables, e and /3, 


as in the case of the barotropic fluid of the non-relativistic fluid. 

When the fluid has a conserved charge, we need to include it as another hydrodynamic variable. 
The conserved current for this charge is a vector and the most general form is expressed as 


ip = mdh 


(14) 


where n is a new proper scalar density. In any frame, n /1=0 should be identified with the net charge 
density. The charge conservation is expressed in the form of the continuity equation as is the case of 
the energy-momentum tensor, 

(15) 


d = 0 . 


If there are more than one charge, we should introduce Eq. (14) for each charge density. But for the 


ideal case, the velocity which defines LRF is universal for all of them, so that the continuity equation 


(15) is valid for each charge density. Therefore, for a system with any number of charges, the continuity 


equations plus one equation which specifies the relation among e, P, and {n^} close the system of 
equations. Normally, as is the case of a non-relativistic fluid, we assume that the thermodynamic 
relation is satisfied for any fluid elements in its rest frame. Then P is expressed as a function of e and 
n through EoS. Note that, in the non-relativistic limit, identifying p m = mn , where n is the number 
density of particles of mass m, and taking f3 <C 1, 7 —> 1 + /3 2 /2, 


(e + P) 7 —> PmC 2 {1 + (pint + P) / PmC~ } • 


(16) 


Here an extra 7 -factor for the energy density is necessary to account for the Lorentz contraction effect 
on the density [S3;, El]- We see that Eq. (12) reduces to Eq. ( 8 ]). 

To obtain the relativistic version of Euler equation Eq. (IT), we first write 


u I yd IJ T' w — ——f huPu^d^uP + hd^vP = 0, 


(17) 


where d/dr = u^d^ = 7 d/dt is the covariant proper time derivative of the fluid element, and h = e + P 
is the enthalpy density. Since u u u u = 1, we have Uyd^u 1 ' = 0. Furthermore, if we introduce the 
thermodynamic relations, 

h = Ts + pn, de — Tds + pdn, (18) 

where s and n are the proper entropy and charge densities, and T and p are the temperature and 
chemical potential, respectively, we can re-express Eq. © as 


Tdfj. ( su M ) + pd^ (nw M ) = 0. 


(19) 


Since n is a proper scalar density satisfying Eq. (15), we conclude that d M (sm m ) = 0, that is, the entropy 
is also conserved. This also holds when the chemical potential p = 0. 


Now, the space components of Eq. (13) in three-vector form is expressed as 


d^T 11 = (/iu) + hud + VP = 0, 

ar 


( 20 ) 


where we introduced the three-vector notation T A ' for { T and u for {u*}, with i = 1,2,3. Using the 


continuity equation for the conserved charges, we can cast Eq. (20) into a more familiar form, 


d 

dr 


e + P 


n 


u = 


-VP, 

n 


( 21 ) 
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which is the relativistic Euler equation. For a system where the contribution to the energy density from 
the chemical potential is negligibly small as in the central rapidity region of ultra-relativistic heavy ion 
collisions, we can use the proper entropy density s in stead of n [95]. 

In the application to the physics of ^/snn > 100 GeV in RHIC and LHC, we often ignore the 
existence of conserved charges and discuss only the behaviors of the energy-momentum tensor when 
we investigate the dynamics of the produced matter at the central rapidity domain. However, in the 
forward rapidity domain, or in the future planned experiments such as BES/BNL, CBM/FAIR and 
NICA/JINR where effects of the large stopping power of nucleus are expected, we cannot ignore the 
contribution from conserved charges in analyses. 

For the sake of the later convenience, we discuss the property of the fluid velocity vP. By using Eq. 
(191), we find 


T^Uy = £U^ 


( 22 ) 


where r is an eigenvector of T and its eigenvalue is e. Physically, this can be interpreted as u M 
being parallel to the energy flow. From Eq. (14), the charge current is proportional to vP. Thus, the 
fluid velocity defined here is parallel to both of the energy and charge flows. Moreover, because of this 
coincidence, e coincides with the energy density when n agrees with the charge density. This is a benefit 
of the ideal fluid approach because we do not have any ambiguity for the choice of LRF and hence the 
application of EoS, differently from the case of viscous fluids. 


3.3 Relativistic Dissipative Hydrodynamics 

For the ideal fluid, we have considered an idealized situation where T^ u and AT are expressed by only 
three quantities, s, n and u^ (five dynamical parameters) due to the assumption of LTE. When we do 
not have these constraints, the general expressions of T^ u and N M contains nine independent components 
more. We expressed them as 


T 1V = (e + P)vPu v - Per + n^, (23) 

j\T = nr + zT, (24) 

introducing the additional quantities, a symmetric tensor Tl^ u and a four vector zT. In order to avoid 
the double counting from already included in the ideal part, these new quantities should satisfy some 
constraints. Conventionally, to be consistent with the usual hydrodynamic form of the Navier-Stokes- 
Fourier (NSF) equation, U 111 ' is decomposed as 

w v = - A^n + rr + rr + (25) 


where n, IP and 7 T" are new variables to characterize the tensor and Ais the projection operator 
orthogonal to iT, defined by A MV = — vPu v . The vectors and tensors appearing here are defined to 

satisfy the following orthogonal conditions, 

u li r 1 ' = 0, hPu^ = u^Ufj, = 7t m m = 0. (26) 

I 11 the ideal fluid, we can choose the scalar functions e and n so as to agree with the energy density and 
the conserved charge density in LRF, respectively, so that LRF can be defined uniquely. Unfortunately, 
this natural identification is not applicable for the above equations. For example, in LRF of the fluid 
flow, where by definition, vP —y (1,0,0,0), 1\T still has the spatial components if zT is not identically 
null: 

AT —> (n,v x ,v y ,v s ). (27) 

That is, in LRF of r, there is still a flow of the conserved charge. Although n = u^N^ 1 is the charge 
density observed in the LRF of r, it is not the same as the density in the LRF of AT. It is also true 
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for e, if h M is not identically null and we cannot identify e with the net energy density measured in LRF 
of the energy flow. 

However, we still have the freedom to choose u 11 to obtain a physically convenient LRF. For example, 
when we observe the fluid velocity with respect to the energy flow, it is natural to employ Eq. (22) as 
the definition of u^. In this case, we should set = 0 and then e can be identified with the energy 
density of the fluid element in its LRF. However if we do so, then in general we cannot set = 0 
simultaneously, and n in the energy LRF defined above does not coincide with the charge density in its 
own flow (current) in LRF of its own flow. The charge density in its own LRF, say n prope r, is given by 


n 


proper 






(28) 


Physically, this modification comes from the Lorentz contraction of the LRF of the charge current with 
respect to the energy LRF since the charge current is not in equilibrium with the energy flow. This 
identification of the fluid flow as that of the energy is proposed by Landau and Lifshitz [93] and referred 
often to as Landau frame. 

On the other hand, we can set v M = 0, instead. Then N M is parallel to the fluid velocity and hence 
n can be identified with the proper charge density in LRF, while £ is not. This definition of the fluid 
velocity is proposed by Eckart |96j and referred to as Eckart frame (see for another choice, called B 
frame, Ref. [97]). 

From the view point of relativity, the difference among Landau, Eckart, /3, or any other frames 
should not cause any physical effects if the theory is truly relativistically covariant. However, since 
hydrodynamics is an effective theory, some approximations or truncations involving thermodynamics 
are implicitly included as we develop more in detail later in Sec. [6] In particular, the thermodynamical 
relations are only defined in a privileged frame, so that its use in a some particular system may lead to 
physically different approximations. In other words, a general C-G method consistent with relativity to 
extract the fluid behavior from a microscopic dynamics has not yet been established as will be discussed 
later in Sec. IS 

In this review, we mostly adopt the Landau frame of the fluid velocity unless mentioned otherwise. 


3.4 Relativistic Navier-Stokes-Fourier model 


In non-relativistic hydrodynamics, an unknown tensor is determined by employing linear irreversible 
thermodynamics [98J. When we apply the same argument to define the corresponding non-relativistic 
quantities to n, n IX1 ' and and write them down in a covariant form, we obtain 


n = -C ax, = 2r]A^d a u^ 


l/ 1 = K 


nT 
£ + P 


a 


where 


2 3 


(29) 


(30) 


is a symmetric direct product of projection operators orthogonal to the velocity field u G The £, 7] and 
k are the coefficients of the bulk viscosity, shear viscosity and charge diffusion, respectively. The above 
expressions are obtained from Ref. [93] . 

This is the relativistic generalization of the NSF theory. In fact, we can derive the NSF theory by 
taking the non-relativistic limit of this model. However, it is known that this model allows a propagation 
which is larger than the speed of light and is essentially unstable. This is discussed in Sec. |3.6 


Sometimes, the above NSF approach is called the first-order theory, since it corresponds to the 


first order derivative expansion of the energy-momentum tensor (see Sec. 6.1.4) or Boltzmann equation 
as is done in Chapman-Enskog formalism (see 6.1.2). Similarly, the causal relativistic hydrodynam¬ 


ics described in the next section is the second order theory or Israel-Stewart type theory, since the 
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corresponding thermodynamics is extended to account for the second order deviation from the thermo¬ 
dynamic equilibrium. As we will see later, it is not known yet that dissipative hydrodynamics can be 
derived from a microscopic dynamics by some systematic expansion from the equilibrium state of the 
system. To avoid misleading, we will not use the terms such as first order and second order theories. 

3.5 Causal models 


In the relativistic NSF theory, the dissipative terms given in Eq. (29) are treated as dynamically 


dependent variables and this is the origin of the inconsistency with relativistic kinematics. In other 
formulations of dissipative hydrodynamics, the number of independent hydrodynamic variables should 
be increased from that of the ideal case given by £, n and wF The new variables are the bulk viscosity 
II, the shear viscosity and the charge diffusion z/h In contrast to relativistic NSF theory, we need to 
construct the evolution of these variables and so far there is no established theory for these dynamical 
equations. There are several variations (see Sec. [6]). 

As an example, we write down the expressions obtained by Israel and Stewart [99i T001. They 
generalized the traditional argument of the linear irreversible thermodynamics and assumed that the 
entropy (four current) has additional contributions attributed to irreversible currents. Following the 
Israel-Stewart theory, let us consider that the entropy four current S 11 is expressed as [991 HDD] 


S ' 1 = So + - Qr 

where Sq is the entropy four current for the ideal fluid and 

Cr = ( — n 2 + ^-TT^TT^ + -I/X 

2 V r n Ty t k 


(31) 


(32) 


To satisfy the algebraic positivity of the entropy production rate, Td^S^ = cr > 0, we find the 
evolution equations of II, zA and 7as 


n = -cax 


dU 1 

Tn~. - -TuUd/jU 1 

dr 2 


TA-n— ( Kl 

2 dr \T( 


zA = 






T]d a u 0 


dVv 1 rj u 

T K ~ - 

dr 2 M 


Tk d_ ( r K 
"dr \Tk 


dr, 


T-k 


a/3 


dr 


1 Q A _ Tjl d_ 

aftOxll ^a/3 i 

A A ar 


r n 

Trj 


(33) 


where d/dr = u^d^. As seen from the above equations, the time derivatives of the newly introduced 
variables are contained in a complicated manner, forming a coupled dynamical system together with 
the original hydrodynamic variables. In addition, besides r), ( and k which have already appeared in 
the relativistic NSF model, other transport coefficients , rn, r K and r n appear in this model. These 
parameters are called relaxation times which originated from relatively rapid dynamics compared with 
that of the relativistic Euler equation, and they play crucial roles to cure the problem of internal 
inconsistencies in the relativistic NSF theory. 

To see this, let us consider an irreversible process of a charged density n satisfying the continuity 
equation, <9tn(x, t ) + V • J(x, t) = 0. Phenomenologically, the irreversible part of the current J is known 
to obey Fick’s law, 

J(x,f) = -PVn(x,f). (34) 

When we consider the retardation effects which are expected from a microscopic dynamics, such a 
current has a structure given by 


J(x,t)= / dsG(t — s)Vn(x, s), 


(35) 
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where G(t) is the memory function which is expressed as the time correlation function of microscopic 
currents. When the memory of microscopic correlations vanishes rapidly, G(t) approximately behaves 
as a Dirac delta function and then we recover Fick’s law. However, when there is no clear separation 
between microscopic and macroscopic time scales, we should not ignore completely the memory effects. 
For this purpose, let us assume that G(t) is approximately expressed as 


G(t) = —e~ t/TD , 
td 


(36) 


where we introduced the relaxation time To, which represents the time scale of the memory effects. 
Substituting this into Eq. (35), we have 


J(x, t) = —DVn(jx.,t) — Tp<9 t J(x, t). 


(37) 


This equation is called the Maxwell-Cattaneo-Vernotte (MCV) equation. If we take the vanishing limit 
of To, Fick’s law is reproduced (see Refs. |101l 1102] for more details). Comparing the above equation 
with equations (33), one can see that rn, r n and r K should be interpreted as relaxation times of bulk 
viscosity, shear viscosity and charge diffusion, respectively. The last term corresponds to the second 
terms of the right-hand sides of Eq. (33). It is also possible to understand the role of the third terms in 
Eq. (33) qualitatively considering that these irreversible currents are densities of corresponding extensive 
quantities (see Ref. [ 103] ). For a more precise derivation of the MCV equation from microscopic theory, 
see Sec. 16.1.31 

Note that in this model, it becomes clear that the values of relaxation times are intimately related 
to suppress superluminal propagation modes. From a linear analysis, we can calculate the propagation 
speeds and in order to maintain relativistic causality, the following relations should be satisfied [ 104 ,. 105 ] 


c 


r n (£ + P ) 


< l-d 


V d 23 

Me + P)- 4 ( s) ’ 


(38) 


where c s is the sound velocity. As we will see later, if the causality is not satisfied the hydrodynamic 
mode becomes unstable, demonstrating that the construction of relativistic fluid dynamics requires a 
special care to unify thermodynamics and the theory of relativity. The former needs a finite space-time 
domain in which a large number of internal degrees of freedom is accommodated, whereas the latter 
demands the true microscopic locality. A similar care should be taken when we deal with a local gauge 
invariance with hydrodynamics. 

The importance of the causal models of dissipative hydrodynamics in relativistic heavy ion collisions 
were first pointed out by Refs. [10614108] . Throughout this review, we use dissipative hydrodynamics 
instead of viscous hydrodynamics, which is commonly used in the literature, since we consider not only 
viscosity but also diffusion (heat conduction). 


3.6 Causality and instability 


The instability of relativistic hydrodynamics has been investigated by various authors. For example, 
Hiscock and Lindblam discussed the stability of causal and relativistic NSF models and concluded that 
the relativistic NSF models both in Eckart and Landau frames are unstable for a linear perturbation 
around a hydrostatic state eds mg. As shown in the previous section, the physical reason for this is 
that the relativistic extension of the NSF theory treats the thermodynamic reaction as a kind of action 
at a distance within the C-G scale. Thus, such a theory necessarily generates acausal modes. 

In Refs. Kill 1105] , it was pointed out that there seems to exist a relation between causality and 
stability: if an acausal mode is contained, such a theory becomes unstable at any cost. For example, it 


was shown that, as far as the causality requirement, Eq. (38) is satisfied in causal models, the stability 
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of linear perturbation is guaranteed in any Lorentz boosted frame. On the other hand, if they are 
violated, the nature of stability becomes Lorentz non-invariant, and unstable mode emerges when the 
hydrostatic background is boosted. That is, the relativistic NSF theory contains acausal modes and is 
inconsistent with the relativistic kinematics. The table [l] shows the stability of various models. 


Table 1: Stability of various models. 


hydrostatic state 
moving frame 


Relativistic 
NSF model 
unstable 
unstable 


causal model 


without Eq. (38) 


stable 

unstable 


causal mode 


with Eq. (38) 


stable 

stable 


The stabilities around Bjorken’s scaling solution are discussed for the relativistic NSF theory m 
and causal model m- Both can be unstable under certain conditions. Moreover, there are various 
works where the relation between heat conduction and stability is discussed. Contrary to the above 
conjecture, there are various proposals to develop stable relativistic NSF theory maum- However, 
the consistency of stability for the Lorentz boost is not studied in these works. 


4 Application to Relativistic Heavy Ion Collisions 

The hydrodynamics described in the previous section as it is cannot be applied immediately to the 
process of relativistic heavy ion collisions. To perform the hydrodynamic approach, we have to specify 
the initial condition (IC), the equation of state (EoS), the transport coefficients, and the freeze-out (FO) 
procedure to connect the hydrodynamical variables to the observed final state particles. Furthermore, 
it is crucial to extract appropriate observables which clearly reflect the hydrodynamic scenario from the 
thousands of produced particles. In below, we describe the necessary ingredients for the application of 
hydrodynamic approach to relativistic heavy ion collisions. 

4.1 Initial Condition 

As shown in the previous section, hydrodynamics is a system of classical held equations with the use 
of thermodynamic relations for each fluid element in its LRF. Initial conditions (IC) for the system are 
the distributions of energy, charge densities, and the velocity fields at an initial time to appropriately 
chosen. For this, we need to know the transition process where the quantum state of two colliding objects 
(nucleus) is converted into the macroscopic matter that defines an initial distribution of hydrodynamic 
variables. This is a hard task since we are not able to perform real ab initio description of collisional 
processes in QCD without introducing models. Dynamics at early stages in the collision is exactly what 
we wish to study with the help of the hydrodynamic approach. As mentioned before, for heavy ion 
collisions, the centrality (impact parameter) is relatively well-estimated on an event-by-event (EbyE) 
basis. Although the centrality is a crucial factor to determine the overlap geometry of the nuclear 
collisions, this is not enough to determine more detailed features of initial conditions. What we have 
to do is to introduce a model which describes the initial collisional stage reflecting the event geometry. 
From this we construct the initial condition of the hydrodynamic equations and then compare the final 
observables in hydrodynamic calculations with the experimental data as function of collision centrality 
and other final state kinematic parameter dependences. 
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Just to understand the qualitative feature of the hadronic collisions at ultra-relativistic energies, 
note that the incident hadrons are predominantly emitted in the forward direction, with high longitu¬ 
dinal momenta (leading particles). The inelasticity distribution P(xf) of the very high energy proton- 
proton (p-p) collision is approximately constant where xf is the Feynmann variable, P(xf) ~ const. 
having relatively small (~ 300 MeV) transverse momentum. This implies that in terms of rapidity 
y = (1/2) In (E + p z ) / (E - p z ), we get P(y) rs-/ coshy, so that the incident protons are emitted pre¬ 
dominantly in the large rapidity domain. Nevertheless, since the average inelasticity is on the order 
of 1/2 ~ 1/5, a large amount of the incident energy is lost and the corresponding energy is used to 
create the particles. The rapidity distribution of these particles are symmetric at the central rapidity 
because most of them are created from the excitation of the vacuum. Therefore, in ultra-relativistic 
p-p collisions, the predominant contribution to the rapidity distribution of charged particles at the 
mid-rapidity comes from the produced particles. See Ref. JB5J for details. 

For A-A collisions, due to the multiple collisions inside the nucleus, the stopping power increases for 
heavier systems, and the baryon number density also tends to immigrate to the smaller rapidity domain, 
especially for lower energies. In the early years hydrodynamic simulations of relativistic heavy ion 
collisions at energies below 1 GeV per nucleon were mainly devoted to see how the nuclear compression 
effects reflect in the observables such as the collective direct flow [63) 114. 115] and full (3+l)D relativistic 
hydrodynamic codes were developed in these works. However, in ultra-relativistic collisions, the baryon 
numbers carried by the colliding nuclei tend to be distributed in the larger rapidity domain while the 
produced matter dominates around the central rapidity. In 1983, J.D. Bjorken m made an estimate 
of the energy density of the produced matter in the central rapidity region of A-A collisions based 
on the boost invariant solution of the (1+1)D longitudinal hydrodynamics. Many of early works on 
hydrodynamic properties from the produced QGP used this kind of initial condition [117J. Some (3+l)D 
codes for such purposes also have been developed (see Refs. [118, 1119] ). 

Later, phenomenological or theoretical studies on hadronic collisions advanced, and several versions 
of more realistic models became available. However, in the early studies of the data from RHIC, 
which started in 2001, most of hydrodynamic calculations were based on smooth energy and baryonic 
distributions (smooth IC). Furthermore, many calculations were concentrated on the analysis of the 
central rapidity region of RHIC data so that (2+l)D calculation has mainly been applied to study 
the transverse dynamics. For these investigations, the smooth IC’s are obtained basically from the 
geometry of the overlapping area of two colliding nuclei. For comprehensive reviews, see Refs. mm 
and references therein. 

More elaborated estimates can be obtained by other microscopic methods based on nucleon-nucleon 
event generators such as HIJING [ 321 ], PYTHIA [122] , NEXUS [ 323 ], UrQMD [ 123 ], EPOS [T25] . 
AMPT [ 126] etc. See Table 3 of Ref. [15]. An essentially different approach is to estimate the QCD color 
field density from the QCD vacuum structure in the high energy limit which is called the Color Glass 
Condensate (CGC) [ 127 , [ 128 ] . In the CGC model for hadronic collisions, the two incident hadrons are 
described as the CGC sheets due to the Lorentz contraction. While they pass through, their interactions 
create the coherent classical gluon field (glasma) which in turn transforms into QGP. Using these models, 
the initial energy-momentum tensor are estimated from the produced gluon number density. Another 
form to introduce the gluon saturation with the mini-jet contribution from nucleon-nucleon collision 
has been developed in the so-called EKRT model. See the recent results in Refs. [ 122 , 150 ] ■ 

Event generator simulations exhibit usually very bumpy structures in transverse energy density 
distribution. Implications of such IC in observables were pointed out by Ref. urn The first use of an 
event generator for the calculation of the initial condition for (3+l)D hydrodynamic simulations was 
performed in 2001 [152 11133] using the Nexus model following the suggestion by K. Werner [12311134] , 
There, the importance of EbyE fluctuations and granularities in IC are emphasized [ 1551 1351 1136] . 
The EbyE fluctuations in IC became the central issue after G. Rolland [137] pointed out that such 
fluctuations in IC are crucial for the analysis of higher order harmonics of collective flow parameters 
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(see discussion in below). 

Presently there are several approaches to generate IC. One is the nucleon-nucleon event generator 
based Monte-Carlo approach which generally called as Monte-Carlo Glauber model (MC-G). There are 
several versions depending on the event generators employed. Instead of using the Glauber approach, 
it is also possible to utilize the microscopic transport models in parton/string level, for example, the 
UrQMD model [ 1241 [1381. I59 J and the PHSD model | H40H143j . Another line is to use the CGC-glasma 
picture, such as MC-KLN ]144H146| and IP-Sat/Glasma |147f4l49j approaches. 

There exists a conceptual difference between the CGC-glasma approach and the other event gener¬ 
ator type approaches, which leads to very different pictures for the initial state for the hydrodynamics. 
In the CGC-glasma picture, the coherent classical gluon held develops longitudinally between the two 
CGC sheets as a consequence of color charge exchange when the two sheets collide. Such a state is 
composed of color flux tubes, which eventually decay into gluons and thermalize, forming the QGP. In 
this picture, the flux tubes are boost-invariant, and particles are generated by the instability of such co¬ 
herent classical color fields. The momentum distribution of produced partons then becomes anisotropic, 
having mainly transverse momenta. On the other hand, the initial condition generated from the par- 
tonic collisions has the large longitudinal momentum as will be shown in Sec. [TJ Such an anisotropy in 
the initial momentum distribution means that the existence of large values of dissipative quantities in 
the corresponding T # “', otherwise the pressure would be isotropic. In most of hydrodynamic simulations 
the initial dissipative quantities are set to null. This constitutes a fundamental problem with respect 
to the initial states of the QGP. 

Furthermore, depending on the initial condition model, the initial profile of the energy density dis¬ 
tribution changes appreciably. Dumitru and Nara pointed out that fluctuations of multiplicities of 
produced particles in p-p collisions within the CGC approach substantially increase the initial fluctua¬ 
tion dominated moments, for example 63 m- 


4.2 Equation of State and Transport Coefficients 

Another fundamental input to execute hydrodynamic calculations are a set of all the thermodynamic 
relations (which in general we refer to as equation of state (EoS)) and the transport coefficients which 
are introduced in Sec. 13.51 

The independent parameters which specify EoS are T and /i#- As mentioned previously, in the 
central rapidity region of ultra-relativistic heavy ion collisions the predominant amount of deposited 
energy goes to the matter created from the QCD vacuum, so that particles and anti-particles are equally 
populated and hence /is is very small. For null baryonic chemical potential, the results of Lattice QCD 


(1QCD) calculations are considered to be well-established [[15111152] as mentioned in Sec. 2.2, even at a 
quantitative level, since some differences among the calculations from different groups, which has been 
existing for a decade, is now tending to converge (see Refs. [771 EDS]). For the temperature well below 
the pion mass, the 1QCD results have been considered less reliable so that the hadronic resonance gas 
EoS is often employed [154]. The 1QCD calculations for finite baryonic chemical potential are non¬ 
trivial due to the so-called sign problem of the fermionic determinant which appears in the partition 
function, but several 1QCD calculations for hnite [ib values have been studied [176] • See more details 
in Refs. [77, 1T55]. A reliable EoS for hnite /z b is crucial to perform a realistic (3+l)D hydrodynamic 
calculations for the heavy ion collisional events to be investigated in BES/CBM/NICA programs. In 
such a hnite (13 domain, effective theories for QCD are often used to obtain EoS [21 EL | 8 l] but a full 
hydrodynamic calculation has not yet been implemented with such an EoS. 

Although the importance of viscosity in the application of hydrodynamics in relativistic heavy ion 
collisions has been pointed out at the early times [ 156] and the necessity of using the causal dissipa¬ 
tive hydrodynamics [106H108] . initially most hydrodynamic analyses of how data at RHIC was with 
ideal hydrodynamics, mainly due to the complexities in the implementation of relativistic dissipative 
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hydrodynamics. Considering the explosive nature of the hot expanding matter in very short time scale, 
it would be rather surprising if LTE is attained in such a system. Nevertheless, the behavior of the 
collective elliptic flow parameter v 2 and other bulk observables (see the subsection below) are quite well 
described by ideal hydrodynamics. The large value of v 2 can be interpreted as the large pressure gra¬ 
dient at the initial stage of the hot fluid expansion, leading to the formation of a new state of matter, 
QGP. Discrepancies between ideal hydrodynamic calculations and v 2 as a function of centrality and 
transverse momentum were expected to be attributed to the effects of viscosity, especially those from 
the hadronic phase where the mean-free path is expected to be much larger than that in QGP. There¬ 
fore in ultra-relativistic heavy ion collisions, the so-called Little Bang PUj picture has been accepted, 
where a hot and dense (almost perfect) fluid of strong interacting plasma of quarks and gluons (sQGP) 
is created together with a rather thin mantle of viscous hadronic matter (core-corona), and eventually 
turns into the bundle of free-streaming hadrons (see Sec. 4.3). 


After the conjecture of the existence of a lower bound called KSS bound of the shear viscosity, 
77 /s > 1 / 47 T where s is the entropy density based on the N=4 Supersymmetric Yang-Mills theory [157 ] 
in 2005, it provoked a great interest in determining quantitatively the coefficient of the shear viscosity 
of sQGP. This is because the success hitherto obtained by the ideal hydrodynamic descriptions of many 
observed collective variables suggests that the shear viscosity should be very close to this conjectured 
lower bound. For a viscous case, we need to specify the transport coefficients such as 77 and (, together 
with the relaxation times, and tjj. Generally they depend on thermodynamical quantities, but in 
the most of calculations, such dependencies are simplified or even taken as constant for 77 /s and r n , and 
bulk viscosity is simply omitted because it is expected to be small. The analysis of the collective flow 
parameters using temperature dependent transport coefficients was introduced by Ref. [158]. For the 
recent studies, see Refs. [ 1301 1159H161| . In addition, the bulk viscosity is expected to be fairly large 
near the critical point, and for a quantitative analysis of experimental data, its effect should not be 
neglected [159 ] [16214168] . 


4.3 Freeze-out and Final State Interactions 

As already mentioned in the Introduction, in relativistic heavy ion collisions it is dificult to observe 
directly the time evolution of flow profiles of the matter produced. Like most explosive evolutions, 
the hydrodynamic description of hot and dense matter eventually fails when the local density and 
temperature become low enough to the point where the continuum description loses its meaning. In 
this stage, the constituent particles gradually decouple from each other and finally turn into a bundle of 
free-streaming particles which are captured in the detectors. We refer to this transition process as the 
freeze-out (FO) process. More precisely, this is called kinetic freezeout in contrast to chemical freeze-out 
for which hadronic chemical composition ceases to change due to the suppression of inelastic channels 
as the temperature and density decrease. Temperature and chemical potential determined from thermal 
models in Sec. prefer to this chemical freezeout. The FO process is a dynamical process and a precise 
description requires an elaborate transport approach as is discussed later. 

In many simplified hydrodynamic calculations, we assume the so-called sudden FO (sFO), where 
the particle spectra is calculated directly from the hydrodynamic variables at the kinetic FO surface. 
However, if we do not employ the chemical FO which occurs earlier than the kinetic FO, abundances 
of particles associated with the pair production mechanism, such as anti-proton, would be totally 
underestimated. The effect of chemical FO is incorporated simply recording the values of temperature 
and chemical potential for each hadronic species at the chemical FO and their abundances are adjusted 
by normalization. However, even in such a simplified picture, changes in the EoS after the chemical FO 
should also be accounted for by treating separately chemically frozen-out particles m- 

Theoretically speaking, the description of the transition from hydrodynamic variables to a bundle of 
free-streaming independent particles should necessarily be of a probabilistic nature, since the former is 
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based on macroscopic distribution functions and the precise microscopic information which is required 
for the later is smeared out from the beginning. In practice, what we can do is at best to calculate the 
final state single particle spectra d 3 a/dff for each observable particle from the hydrodynamic calculation. 

In order to introduce a realistic FO in a hydrodynamic calculation we first have to know when FO 
starts, i.e. determine what is the physical criteria to switch the hydrodynamic evolution to some non¬ 
equilibrium transport scheme. Second, we need to map the hydrodynamic quantities (fluid) to a state 
of interacting many particle system. Third, a dynamical description (transport equations) of this many 
particle system should be implemented till all the particles become free-streaming. After this stage, we 
still consider processes of decays of unstable resonances down to the observable mesons and baryons. 

Each of these problems is non-trivial, which constitutes formidable work. The simplest proposal is 
the previously mentioned sFO from Cooper and Frye [33] . which neglects the third step above. In the 
sFO, we first establish the so-called freeze-out surface (FOS), E, using some physical criteria. Many 
people uses a constant FO temperate T FO for simplicity but other ways have also been proposed H2B- 
1176] . In Fig. [3j a typical time evolution of FOS in the transverse direction at central rapidity for a 
smooth initial condition is shown. 



r [fm] 


Figure 3: (color online) Example of isothermal hypersurfaces. 


Figur^] taken from Ref. m- 


Cooper and Frye proposed to calculate the freeze-out hadronic spectra as 

d 3 N T 

= / s d(Tfl pl ‘^ l) ’ ( 39 ) 

where (x,p) is the distribution function of the specie i at FOS. In the case of an ideal hadronic 
resonance gas (HRG), the distribution function, /, is given by the Maxwell-Jiittner-type distribution 
function (see the later section on the Boltzmann Equation) at thermal equilibrium, 


f(x,p) ->• f TE (x,p) 


9 _ 1 

(2nh) 3 e( u -P~ri/ T ± 1’ 


(40) 


where g is the statistical factor of the particle and the sign + and — is for bosons and fermions, 
respectively. Differences in hadronic species enter through g,p and its mass [12] , This expression is a 
natural choice for ideal hydrodynamics where LTE is satisfied. However, in the presence of dissipative 
effects, the distribution / is not given by Eq. (40) but there exists a correction, 


/ {x,p) = /te (x,p) + A/ (x,p ). 


(41) 


“Reprinted with permission from H. Niemi and G. Denicol, How large is the Knudsen number reached in fluid dynamical 
simulations of ultra-relativistic heavy ion collisions?, arXiv:1404.7327 (2014). 


20 





















When the dissipative corrections are small, A/ can be related to rP 1 ' as well as the thermodynamical 
quantities. We will discuss more in detail in Sec. 6 . 1.2 together with the derivation of dissipative 
hydrodynamics from the Boltzmann equation. 

Unfortunately, using the Cooper-Frye formula (CF) one cannot directly calculate the identified 
particle spectrum, since the FO temperature for chemical FO (that is used for the thermal model for 
particle ratios) is higher than hydrodynamic FO (kinetic freeze-out) [ TT81 Il79lj . However, the dominant 
part of the produced particles is pions and chemical non-equilibrium effects from heavier baryons are 
considered to be small. Thus, the modification of the EoS due to the chemical FO is often neglected. 
Even in such a simplihed treatment, many hadronic resonances need to be taken into account to calculate 
the correct multiplicities of observed identified particles through their decays down to the stable hadrons 
at the kinetic FO temperature, which is usually taken around Tpo> 120 MeV. Several open codes which 
describe the thermal HRG are available and can be used to obtain the final stable hadrons from the 
sudden FO hadronic states [87, 88| Il80fil82| . For more discussion on this point, see Chap. 26 of Ref. 


Another shortcomings of CF is the so-called negative contribution problem. As seen from Fig. [3] 
FOS has two distinct regions: one in which the normal vector to £ is time-like (that is, the surface 
itself is a space-like 3D volume), and the other space-like (the surface is consisted of 2D space surface 
expanding in time). Particle emissions in the domain where the normal is time-like (volume emission) 
has no conceptual problem in CF. On the other hand, in the domain where the normal is space-like 
(surface emission), a portion of particles in the distribution / belongs to the fluid, and should not be 
counted as emitted particles from the fluid. However, if we only count the particles which are going out 
of the fluid surface, then the total energy of the FO particles does not coincide exactly with the total 
energy of the fluid [BPIj . 

In order to improve the problems of sFO, the so-called hybrid models have been developed where 
the final stage of hydrodynamics is connected to the event generator based hadronic cascade codes 
or Boltzmann type equations [2D]. Importance of hybrid approach has been demonstrated by various 
groups. One of the advantages of these approaches is that the chemical FO processes is naturally 
incorporated in the scheme. More detailed description of the hybrid model, see the recent reviews by 
Hirano et al. [15] and by H. Petersen m■ 

In most hybrid approaches the transition from the hydrodynamic description to a hadronic transport 
system is done using CF. Therefore, the problem of CF mentioned above still persists, since the switching 
surface has the similar characteristics as the final FOS. In fact, a small correction originated from the 
negative contribution may be masked by the fluctuations and uncertainties in the proceeding hadronic 
cascade processes, but this may not be the case for other observables such as HBT radius (see the next 
section) especially for Kaons, and particle ratios in smaller systems such as p-p and p-A. In addition, 
the negative contribution problem is also a question of theoretical consistency and needs a solution. 
Several approaches have been proposed and discussed in |183fH8811188] . The method proposed in Ref. 
[ 188] has been applied to describe these observables and shown to be effective [18911194] . One important 
element of such an approach is that when we switch from the hydrodynamic phase to the hadronic phase, 
the EoS of 1QCD and that of the hadron resonance gas should be consistent otherwise the conservation 
of energy and momentum is compromised (see Ref. [ 195] ). It is also pointed out that the hadronic mass 
spectrum affects the final state collective flow parameters [196]. 


4.4 Physical Observables 

The signals related to the hydrodynamic evolution in the experimental measurements are those associ¬ 
ated to the collective flow profile. In hydrodynamic regime, if the initial condition has spatial anisotropy, 
the corresponding pressure gradient will generate the acceleration of the fluid elements, which in turn 
will affect the angular distribution of the final particles. On the other extreme, if the produced particles 
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come from just a superposition of independent nucleon-nucleon collisions and there is no interaction 
among the produced particles, then the final angular distribution is just a sum of single nucleon distri¬ 
bution and the initial collision geometry will not be reflected. In this sense, an evidence of a systematic 
dependence of final particle distribution on the initial anisotropy indicates strongly the existence of 
some collective dynamics. However, anisotropy in the initial condition is not directly observable. Even 
so, if it is related to that created by the non-central collision geometry, it should be strongly correlated 
to the multiplicity of the produced particles (centrality). 

To more quantitative studies, the Fourier expansion of the transverse angular distribution is most 
commonly applied prHigg] . For very high multiplicity events, we can introduce the differential particle 
distribution for each collisional event as 


d 3 N 


d 2 N 


PxdPxdyd(p 2n pxdpxdy 


1 + 2 J+„cosn (<f> - i/>rp) 


(42) 


n =1 


where ip up is the angle of the reaction plane (the plane defined by the incident direction and the impact 
parameter) of the event with respect to the observational frame, i.e. the azimuthal angle of the impact 
parameter b of the collision in the observational frame. In the above equation, it is assumed an ideal 
situation where the initial condition has the smooth distribution and symmetric with respect to the 
reaction plane (RP). The Fourier coefficients, v n are called collective flow parameters and are functions 
of px and y. Sometimes iq is called as directed flow, v 2 as elliptic flow, and v 3 as triangular flow. 
For smooth initial conditions all the odd flow parameters should vanish at the central rapidity region 
(y ~ 0), and for the central collision (b ~ 0) all the flow parameters should vanish. 

Before 2010, most hydrodynamic calculations of elliptic flow for the RHIC experiments were done 
using smooth initial conditions, except for those mentioned in Sec. 4.1 In a theoretical calculation 


we know the event plane a priori, but experimentally we are not able to control the reaction plane 
in each collisions. Furthermore, each event does not have smooth and symmetric behavior, so that 
the comparison of theoretical values to experimental ones requires certain care. Several methods to 
determine experimentally v 2 as function of centrality and pt have been developed nanssHM]. 


When we need to consider EbyE fluctuations, Eq. (42) is not applicable because there is no symmetry 


on EbyE basis. The most general form of angular distribution can be parameterized as 


d 3 N 

prdprdydqb 


1 d 2 N 
27r pxdpxdy 


OO 


1 + 2 v n cos n (0 — ip. 

n=l 


EP\ 


(43) 


The angles should be determined for each value of px and y so as to make v n maximum. However, 
due to the insufficient number of particles in each bin of px and y, the direct determination of these 
parameters for one single event is not possible. Thus, instead of determining directly these parame¬ 
ters, particle correlation measurements are used to obtain the moment distributions of collective flow 
parameters. For more detailed discussion on the parameters v n and their experimental determination, 
see Refs. [TUI H3 1199] . The event angles are also not measurable directly for each single event but 
correlations among flow parameters can be determined. Such flow-plane correlations are argued as a 
good signal of hydrodynamic evolution [202] . 

Describing the characteristics of collective dynamics generated from the initial state profile is the 
central issue of the hydrodynamic approach, however, the other basic bulk behaviors of the data such as 
identified single particle spectra and the size of emission domain should be reproduced simultaneously. 
The size of emission domain is measured from the so-called HBT interferometry of two identical particles. 
Use of the interferometry of two photons emitted from a stellar surface to measure the size of stars was 
proposed by Hambury Brown and Twiss in 1954 |203] and this method was first applied in multiple 
pion production processes in 1960 by G. Goldhaber et al. [ 204j . Since then HBT interferometry became 
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one of the important methods in hadronic collisions and in relativistic heavy ion collisions in particular. 
For comprehensive reviews, see Refs. |205H212] and the references therein. The application of HBT 
interferometry for pions produced in relativistic heavy ion collisions determines three parameters, i?i ong , 
.Rside and R out which characterize the shape of the emission area of the pions as functions of the average 
momentum of the pion. Here, i?i ong represents the longitudinal radius, R S id e the side direction radius, 
and R 0 ut refers to the radius measuring the optical depth of the emission area. In contrast to the stellar 
case, the emission surface is dynamical and there exist other mechanisms which create two-particle 
correlations, thus, special care should be taken to define these parameters and compare the theoretical 
values to the data [213] . 


5 Overview of the Present Status 

A large number of experimental and theoretical studies on collective flow phenomena and extraction of 
QGP properties have been done in the last 40 years [6TH01 IT2l T3j 1214] . In particular, the concept of 
the formation of sQGP and its subsequent hydrodynamic expansion (Little Bang picture) of relativistic 
heavy ion collisions has been established in the last decade as was announced in the BNL news in 
2005 [215] . It was argued by Hirano and Gyulassy [ 216 ] that the observed large elliptic flow v 2 and 
consequently almost perfect fluidity of the newly created matter can be understood as a sudden increase 
of the entropy in the QGP phase contrasted to the hadronic gas, thus it is considered as a signal of 
deconfinement. 

As already pointed out, in the hydrodynamic picture the observed flow pattern should reflect the 
properties of the initial condition. Therefore, the hydrodynamic description allows for a more detailed 
study of the nature of the initial conditions using now measurable EbyE fluctuations in flow pattern 
and their correlations |1601[202ll217f!2l9] . It has been shown that the behavior of these higher harmonic 
coefficients as functions of particle mass, centrality and transverse momentum are well reproduced by 
hydrodynamics. The hydrodynamic modeling to describe the experimental flow parameters on EbyE 
basis is considered to be a unique tool to study the initial QCD dynamics of nucleus-nucleus (A- 
A) collision at ultra-relativistic energies. Additionally, there are attempts to determine the transport 
properties of QGP from hydrodynamic analysis. As mentioned in the Introduction, there are already 
many recent review articles on these subject written by authors who themselves played important roles 
in these Endings [151 1T7H2U1 [22011221] . Therefore, we refer readers to these references for the development 
of the model and its understandings. 

5.1 News from recent LHC and RHIC Data 

RHIC has been exploring nuclear matter at extreme conditions over more than a decade now and a huge 
amount of data have been collected for collisions of light and heavy nuclei at several energies, ranging 
from below 10 GeV up to the maximum of 200 GeV per nucleon pair in the center of mass frame. 
While the understanding of the experimental data has become deeper and many questions concerning 
the collective aspects of the evolution of the created matter in such collisions have been answered, many 
other questions have also emerged. Subsequent data, in particular from the recent results from LHC and 
the updated RHIC experiments seem to consolidate more firmly the scenario of Little-Bang. Similarly 
to the analysis of correlation measurements in modern observational cosmology [222] , the flow pattern 
characterized in terms of higher Fourier components became possible on a real EbyE basis due to large 
multiplicities in LHC data. 

With the start of the heavy ion program at the LHC in 2010, very precise new data has been 
measured for collision energies of one order of magnitude higher than the maximum achieved at RHIC. 
The rich collection of observables measured by the large and complex experiments running at RHIC 


23 
















and LHC gives us an unique opportunity to study in detail the aspects of the hot and dense matter 
created in heavy ion collisions. In the following we highlight the latest experimental results as well as 
the relevant constraints they are setting for the theoretical developments. 


5.1.1 Anisotropic Flow 


One of the most important signals of collective behavior during the system evolution in heavy ion colli¬ 
sions is the azimuthal anisotropy of the produced particles. Initially, all the analyses were concentrated 
on the elliptic flow parameter, the second harmonic of the Fourier decomposition (see equation 42). 


Since non-central collisions form the so-called almond shape in the overlap region of the two incident 
nuclei, it is reasonable to expect the observation of the elliptic anisotropy in the final state if collectivity 
is present during the evolution [TJ [9j [TO], 1197 . 199]. 

The measurements of the tq coefficient helped to improve our understanding for the collective aspects 
in heavy ion collisions. From the experimental side special efforts have been made to develop techniques 
to compute flow harmonics from the final particles azimuthal distribution and to extract the possible 
influences from non flow and fluctuations (198112011I223H227] . In the last few years higher order flow 
harmonics have attracted special attention as important observables of the structures of the initial 
condition [ 1371122811232] . There has been a wealth of recent results of v n reported by RHIC and LHC 
experiments seen in (2271122911232H23411234H249] . 

Regarding the experimental methods used to extract the Fourier coefficients in anisotropic flow 
analysis, the ATLAS collaboration has recently presented experimental measurements of the EbyE 
v n distributions (see Ref. [244] for details). Such measurements allow to directly extract the mean 
( v n ), the standard deviation cr Vn , and the relative fluctuation cr Vri /(v n ) of the v n distribution, which 
were previously usually only estimated from second and fourth order cumulants assuming cr Vn <C ( v n ) 
[266]. Figure [4] shows the probability density distributions of the flow coefficients obtained from EbyE 
calculations for several centrality windows. The comparison with results obtained with more traditional 



Figure 4: Probability density distributions of event-by-event tq (left panel), tq (central panel) and tq 
(right panel) measured by the ATLAS experiment. Figur^] taken from Ref. [244] . 

methods such as second and fourth-order cumulants (y n { 2}, u n {4}) and event-plane (u n {EP}) present a 
clear systematic ordering between the different techniques, with v n {2} > v n {EP} > v n {EbyE} > v n {4} 

“Reprinted from G. Aad et al ., Measurement of the distributions of event-by-event flow harmonics in lead-lead collisions 
at yA/viv = 2.76 TeV with the ATLAS detector at the LHC , J. of High Energy Phys. 2013 (2013) 183, under the terms 
of the Creative Commons Attribution Noncommercial License. Copyright © 2013, CERN, for the benefit of the ATLAS 
collaboration. 
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for the centrality range measured (cf. Fig. 10 in Ref. [227]), with v n {EbyE} being the mean value 
extract from the EbyE v n distributions. Although non-flow correlations may affect more strongly 
the lower order cumulants (see discussion in | 200 j ). the current experimental techniques are able to 
largely suppress their contribution. Therefore, the ordering observed may be attributed mainly to flow 
fluctuations [1341 11361 1 13911250fl252j . Such effect can be explicitly shown for the second and fourth 
order cumulants, where the finite variance of the v n fluctuations contributes positively to v n {2\ and 
negatively to u n {4} [ 253] . 

Assuming that the development of anisotropic flow is driven by the pressure gradients defined by 
the structures of the initial density profile, the measurement of flow fluctuations may reveal important 
properties of the initial state. However, it is important to remember that the fluctuation spectrum 
extracted from the measurements of flow coefficients can be distorted by nonlinear effects during the 
collective dynamic evolution of the system. The precise connection with the initial fluctuation spectrum 
can be established only when the non-equilibrium dynamics of the system is understood, but such 
measurements will help to constrain the developments of initial condition models [218] . In particular, 
the measurements of flow coefficients using different order cumulants present different sensitivity to 
flow fluctuations and allows for further exploration of the details of the underlying flow distribution. 
Results from PHOBOS [254] and STAR [ 236 ] Collaborations for Au-Au collisions at sJsnn = 200 
GeV are consistent with a Bessel-Gaussian functional form for the u 2 distribution [255]. which implies 
u 2 {4} w 2 { 6 }. Furthermore, very precise measurements performed by the ALICE experiment [256] also 

seem to qualitatively follow a Bessel-Gaussian model of flow fluctuations at the LHC (see the left panel 
in Fig. [5]) . Although, it was found that u 2 {4} is slightly different than u 2 { 6 }, which raised the possibility 
of non-Bessel-Gaussian fluctuations, in agreement with the recently proposed power law distribution 
[257]. The ATLAS Collaboration has also reported measurements of v 2 {6}/u 2 {4} and n 2 {8}/n 2 {4} 
where significant deviations are observed for most central and peripheral cases (see the right panel in 
Fig.g [227]. There have been some extensive investigations on the correlation of the observed flow 
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82 65 52 43 31 17 7 




Figure 5: Measurements of u 2 { 2 }, v 2 {4} and n 2 { 6 } from the ALICEj^] experiment (left panel) and ratios 
of u 2 {6} and u 2 { 8 } to v 2 {4} from the ATLAS^] experiment (right panel) as a function of centrality in 
Pb-Pb collisions at yj snn = 2.76 TeV. Figures taken from Refs. [256] and [227], respectively. 


“Reprinted from B. Abelev et alMultiparticle azimuthal correlations in p-Pb and Pb-Pb collisions at the CERN 
Large Hadron Collider, Phys. Rev. C90 (2014) 054901, doi:10.1103/PhysRevC.90.054901, under the terms of the Creative 
Commons Attribution 3.0 License. Copyright © 2014, CERN, for ALICE Collaboration. 

^Reprinted from G. Aad et al., Measurement of flow harmonics with multi-particle cumulants in Pb+Pb collisions at 
yj Snn = 2.76 = 2.76 TeV with the ATLAS detector, Eur. Phys. J. C74 (2014) 3157, under the terms of the Creative 
Commons Attribution 4.0 License. Copyright © 2014, CERN, for the benefit of the ATLAS collaboration. 
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fluctuations with the initial state profile [SHE, 1258H26U] and, as reported by the ATLAS experiment, 
the eccentricity distributions calculated with Glauber and MC-KLN models [ [261 . 1262] (see Sec. 4.1) 


cannot completely describe the observed data. Another interesting consequence of EbyE fluctuations 
recently pointed out in Refs. [ 2631 !264j is the possible factorization breaking of two-particle correlation 
probability distributions into a product of single-particle distributions. Figure [6] shows preliminary^] 
results of the ratio u 2 {2}/u 2 [2] for differential {jpp dependent) flow coefficients obtained from Pb-Pb 
collisions at y/s^N — 2.76 TeV by the ALICE experiment [265] . The numerator (u 2 {2}) in the above 


ratio is the usual second order cumulant flow coefficient (see for instance Refs. [ 200 :. 201 . 266] ). in which 
only one of the particles is taken from the pr bin of interest, while the denominator (i’ 2 [2]) is computed 
by taking both particles from the pt bin of interest [264]. The factorization breaking (deviations from 
unity) as a function of the transverse momentum is shown to be the result of decoherence in the angular 
correlations induced by initial fluctuations and is predicted even for pure hydrodynamic calculations 
The intensity of the factorization breaking observed in the experimental data is thus of special 


interest in constraining the initial state and transport properties of different model calculations [264]. 
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Figure 6: (color online) Ratio of u 2 {2} over u 2 [2] as a function of pt for three centrality classes mea¬ 
sured by the ALICE experiment for Pb-Pb collisions at ydv)v= 2.76 TeV [265] . Comparison with 
hydrodynamic calculations [ 264 ] are also shown. Figur^] taken from Ref. [265]. 


With the latest experimental results the interpretation and understanding of the enormous amount 
of information contained in particle correlation analysis have become a very challenging task. Therefore, 
newer and more sophisticated techniques need to be developed in order to help us disentangle the true 
meanings of our measurements. In this sense, a very recent work by Bhalerao et al. [267 ] proposed 
the application of the statistical technique of principal component analysis (PCA) to study EbyE 
fluctuations in two-particle correlations (2PC). The method approximates the Fourier coefficients of 
the pair distribution, 14 a, as a sum of flow fluctuation components (modes), 


v,a(pi,P2) 


(44) 


a=l 


where the coefficients 14 a are determined by the statistics of the single-particle complex Fourier coeffi¬ 
cient, T4a(pi,P2) = (V n (pi)V* fa )), and pi is a shorthand notation for the coordinates Pt and 77 [267]. 
In the absence of flow fluctuations the left hand side of Eq. (44) factorizes, i.e. only a = 1 is present, 


thereby recovering the usual anisotropic flow. One of the advantages of this new method is the ability 
to make use of all the information contained in 2PC, therefore, allowing for a more robust observable to 


2 The effects of non-flow contribution on r» 2 {2} and ^[2] are being further investigated. 


“Reprinted with permission from Y. Zhou (for the ALICE Collaboration), Searches for dependent fluctuations of 
flow angle and magnitude in Pb-Pb and p-Pb collisions , Nucl. Phys. A 931 (2014) 949. Copyright © 2014 The Authors. 
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study flow fluctuations. Indeed, it has already been shown that previously unknown subleading modes 
in both rapidity and transverse momentum are found to be present in data measured by the ALICE 
experiment |267l 1268] . which in turn are argued to be responsible for factorization breaking as observed 
in two-particle correlation analysis [2631 126411269] . 


5.1.2 Two-Particle Correlations 


The study of two-particle correlations (2PC) has become one of the main topics of investigation in the 
last years. The general idea consists of constructing a correlation function C(A0, A//) by taking the 
distribution of pairs relative differences A0 and A?/ between a trigger and an associated particle in the 
“same” event and then dividing by the distribution of pairs of particles from “mixed” events, in which 
the trigger is taken from one event and the associated is taken from a similar but different event, which 
thus cancels all the effects due to detector acceptance and efficiency. Given the desired ranges for the 
transverse momentum of the trigger and associated particles, the correlation function can be written 


as, 


C( Act), An) 


Admixed ^ A same ( A0, A'/^) 

Asame A m ixed(A0, A'/^) 


(45) 


As a matter of fact, we recall the observation of structures such as the long-range near-side (A <f> ~ 0) 
ridge and the away-side (A</> ~ n) double-hump in the A g x Acft distributions obtained from Au-Au 
collisions at RHIC [ 27011271 ] . Likewise, recent measurements at LHC also confirm the observation of 
such structures in Pb-Pb collisions [2301 243 . 200 ]. 

The complete understanding of the origin of such structures is still under intense investigation. 
Initially, the proposed interpretations were based on the medium response and conical flow induced by 
jet quenching [ 27211273 ] as well as longitudinal flux tubes in the initial state [ 27311277 ] . However, it was 
later shown that EbyE hydrodynamic simulations with inhomogeneous initial condition could produce 
similar structures [ 278 ] suggesting that the ridge and the double-hump structures could be explained 
by the presence of higher order flow harmonics [ 137 '. 127911281 ] . which appears naturally as the odd 
harmonics contribute constructively at A <f> ~ 0 and destructively at A <f> ~ n. The Fourier coefficients 
in this case are given by, 


14a = (cosnA0) 


Sm=i co 

EL, C(AM 


(46) 


where the correlation function C (A0) is scaled by the respective number of pairs in the mixed and same 
events distributions of the projected Ag slice. 

In the earlier studies of 2PC performed at RHIC for p-p and d-Au collisions, the observed results 
were very similar |282| and the structures observed for Au-Au collisions were not present. Therefore, 
these results were commonly taken as reference, where no QGP medium was formed, to compare with 
Au-Au collisions data. However, recent measurements at the LHC for high multiplicity events in p-Pb 
[28311285] and even in p-p [28611288] collisions have indicated the appearance of the same structures 
usually observed only in A-A collisions. Figure [7] shows a comparison between the two-dimensional 
distributions of charged particle pairs measured by the CMS experiment for low (left panel) and high 
(right panel) multiplicity events [283] , where the ridge-like shape can be clearly observed in the second 
case. 

In addition to the already surprising similarities between the measurements of peripheral Pb-Pb 
collisions and the high multiplicity selection in p-Pb collisions for the 2-dimensional Ag x A <f> correlation 

“Reprinted from S. Chatrchyan et al., Observation of long-range, near-side angular correlations in pPb collisions at 
the LHC , Phys. Lett. B718 (2013) 795, under the terms of the Creative Commons Attribution 4.0 License. Copyright 
©CERN, on behalf of the CMS Collaboration. 
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Figure 7: Arj x A<f> distributions for charged particle pairs measured by the CMS experiment for low 
(left panel) and high (right panel) multiplicity events of p-Pb collisions at = 5.02 TeV. Figure^] 

taken from Ref. |283j . 


function, comparisons of the first few coefficients extracted from Fourier decomposition show remarkable 
consistence in shape between the cases. Figure [8] shows the comparisons between the results obtained 
by the ATLAS experiment for peripheral 55-60% Pb-Pb collisions at ^snn = 2.76 TeV and high 
multiplicity p-Pb collisions at y/smf = 5.02 TeV |289j . The plots on the left show the original values 
extracted from each analysis and the plots on the right show the results extracted from Pb-Pb collisions 
rescaled both vertically, by an empirical factor of 0.66, and horizontally, by a factor of 1.25 as proposed 
in Ref. [ 290] to take into account the difference between the ( pp) values in the two systems. 

Following the observations made by the LHC experiments for the p-p and p-Pb collisions, the 
PHENIX collaboration presented new t>2 measurements for central d-Au collisions at yj snn = 200 GeV 
compatible with the values reported for p-Pb collisions at y/sNN = 5.02 TeV [291] . 

The unavoidable question is whether the dynamics governing the evolution of the system in each 
case, namely A-A, p-A and p-p collisions, is similar or not. If there is collectivity in p-A and p-p colli¬ 
sions, can hydrodynamics still be applied to such small systems? Indeed, elliptic flow measurements of 
identified hadrons in p-Pb collisions by the CMS [292] experiment seem to be consistent with the quark 
number scaling previously observed in A-A collisions [234112351 293 1 294]. Thus, these newly produced 
experimental results have undoubtedly imposed serious challenges to our theoretical understandings. 
The assumptions generally used in the hydrodynamic approach for A-A collisions, namely the thermal- 
ization time, local thermal equilibration and coarse-graining scale, may need to be revised if one wants 
to describe small systems such as p-A and p-p. At the same time, it is a great opportunity to probe 
the edges between the two worlds usually described independently with macroscopic and microscopic 
theories (see interesting discussion in Ref. [ 295] ). 

5.1.3 Direct Photons and Pion Femtoscopy Measurements 

Another important experimental probe of the properties of the hot and dense medium created in rel¬ 
ativistic hadronic collisions is the measurement of direct photons. Since photons are not affected by 
strong interactions they can emerge from the fireball almost undisturbed while carrying information 
about their emission source. 

“Reprinted from G. Aad et ai, Measurements of long-range pseudorapidity correlations and azimuthal har¬ 
monics in yj Snn = 5.02 TeV proton-lead collisions with the ATLAS detector , Phys. Rev. C90 (2014) 044906, 
doi:10.1103/PhysRevC.90.044906, under the terms of the Creative Commons Attribution 3.0 License. Copyright © 
2014 CERN, for the ATLAS Collaboration. 
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Figure 8: v n coefficients as functions of px for Pb-Pb collisions at y/s^N = 2.76 TeV and p-Pb collisions 
at y/sw n = 5.02 TeV. Figur^“] taken from Ref. [289] . 


In heavy ion collisions, direct photons can be classified as prompt photons , which are mainly produced 
by leading order pQCD processes in the very initial stages, and thermal photons , emitted by radiation 
of the medium, both in partonic and subsequent hadronic phases. The thermal photons spectrum 
should therefore depend on the temperature of the respective production source [2961 1297] , Recent 
hydrodynamic modeling suggests the possibility of extracting information of the early hydrodynamic 
evolution of the system (see Ref. [ 298 and references therein). 

Experimentally, the measurement of direct photons is a very complex task that requires careful 
identification and removal of photons from hadron decays in order to isolate the direct contribution 
[292]- Thus, experimental techniques have had to develop over detector upgrades and experiments at 
RHIC and LHC have now been able to present crucial direct photons measurements. 

The PHENIX experiment at RHIC recently observed an exponential excess in direct photon trans¬ 
verse momentum yields for Au-Au collisions compared to p-p collisions [ 299 . 300 ]. The left panel in 
Fig. | shows the combined px spectra for both p-p and Au-Au collisions at y/smf = 200 GeV. The 
excess obtained from the difference between the Au-Au results and the V co n-scaled p-p results is well 
described by an exponential distribution with an average inverse slope of 240 MeV in the px range from 
0.6 to 2 GeV/c [299 ]. In a similar analysis, the ALICE experiment at LHC also observed an excess of 
direct photons in Pb-Pb collisions at yW/v = 2.76 TeV with respect to NLO pQCD predictions scaled 
by iV coll [ 297 ] (see right panel in Fig. [ 9 ]). The inverse slope parameter extracted from the exponential 
fit in the p T range from 0.8 to 2.2 GeV/c is on the order of 300 MeV [297] . 

“Reprinted with permission from A. Adare et al., Centrality dependence of low-momentum, direct-photon production 
in Au+Au collisions at y/s NN = 200 GeV, arXiv: 1405.3940 (2014). 

“Reprinted with permission from M. Wilde (for the ALICE Collaboration), Measurement of Direct Photons in pp and 
Pb-Pb Collisions with ALICE, Nucl. Phys. A904-905 (2013) 573c, Copyright © 2013 CERN. 
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Figure 9: pr spectra of direct photons measured by the PIIENIX^] experiment for Au-Au and p-p 
collisions at n = 200 GeV [ 299 ] (left panel) and results from the ALICEN experiment for Pb-Pb 
collisions at = 2.76 TeV [2H7J (right panel). Figures taken from Refs. [299] and [557], respectively. 


In a hydrodynamic picture, these results can be interpreted as the effective temperature of a ther- 
malized source convoluted over its whole time evolution. In both cases, i.e. Trhic ~ 240 MeV and 
Tlhc rs-/ 300 MeV, the values are well above the phase transition temperature 2% of approximately 170 
MeV predicted by lattice QCD calculations for the transition to the QGP phase [29TJ. Such high values 
of temperature may indicate the production of direct photons at very early times during the system 
evolution [301]. On the other hand, PHENIX and ALICE experiments have also reported measurements 
of significant elliptic flow of direct photons in the low-p^ region for collisions of Au-Au at -\Jsnn — 
200 GeV [1502 ] and Pb-Pb at = 2.76 TeV [501]. respectively. While one would expect a small 

elliptic flow if the direct photons are emitted predominantly at the very early stages, since flow develops 
towards the final stages the observed large values of u 2 also support relevant photon production later 
in the evolution. Therefore, the simultaneous description of the above observations has become truly 
challenging to the current theoretical developments (see Ref. ra). 

Comparisons from predictions of EbyE hydrodynamic calculations [304] as well as from microscopic 
dynamic model simulations using the PHSD model [ 3GB ] have been presented for ALICE and PHENIX 
results, respectively. As shown in these works, the initial condition fluctuation in EbyE hydrodynamic 
calculations seems to play an important role in the final direct photon u 2 [ [304] . Investigations using 
a (3+l)D hydrodynamic model which takes into account viscous corrections suggest that the flow of 
direct photons may be also especially sensitive to the equation of state and the transport properties of 
the system 3~(T6 . Furthermore, it is worth noticing the excellent agreement of the predictions given by 
the PHSD model with respect to the PHENIX measurements (cf. Fig. 11 in Ref. [305] ). Within the 
approach employed in this model, the photons produced during the QGP phase account for slightly less 
than 50% of the observed spectrum and have small u 2 (on the order of 2-3% only), whereas the main 
contribution to the large u 2 comes from hadronic channels. These channels, however, indirectly require 
the strong interactions of the partonic phase to effectively build up the observed asymmetry of the final 
particle momentum distribution. In another recent study employing the hydrodynamic framework [507] 
it was suggested that the precise measurement of direct photons in p-A collisions may furnish relevant 
information on the fluid dynamic profile at early times, in particular, concerning the possible formation 
of a QGP phase. 

Finally, femtoscopy measurements based on quantum-statistical interferometry of identical bosons 
as proposed by Hanbury Brown and Twiss [20311308] . sometimes also called HBT correlations, bring 
additional constraints to the modeling of the dynamic evolution of the system produced in high energy 
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collisions. As discussed previously, the main goal of this experimental technique is to infer about the 
space-time scales of the emitting source [2051 imi 1212] . which allows one to extract information about 
the duration of expansion and the system size at the kinetic FO. The details of the analysis performed 
on the experimental data can be found in the Refs. [3091 - 1312] . Figure 10 shows the dependence of 
the product of the HBT radii (left panel) and the decoupling time 77 (right panel) with the charged 
particle multiplicity for several femtoscopy measurements at different collision energies. The later is 
an estimated quantity obtained by assuming that Ai ong is proportional to duration of the longitudinal 
expansion and, therefore, to the decoupling time of the system (see Ref. [3T2] for more details). The 
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Figure 10: Product of HBT radii extracted from pion correlations at kp = 0.3 GeV/c (left panel) and the 
decoupling time extracted from -Ri on g(^T) (right panel) as functions of the charged particle multiplicity 
for several collision energies. Figure^] taken from Ref. [312] , 


common linear behavior shown in this figure, that extends over such a large range in collision energies 
(also reported in Ref. 310; ). is consistent with a universal condition for pion FO |313| and establishes 
important constraints for model developments. Extended measurements using three pion correlations 
[ 311] and the cumulant expansion [ 314 J have also been recently reported by the ALICE collaboration 
for Pb-Pb collisions at sJsnn = 2.76 TeV as well as for p-Pb at sJsnn = 5.02 TeV and p-p ai^/sNN = 
7 TeV. These improved measurements restrict even further the model developments. In particular, as 
concluded in Ref. [ 314] . the results obtained seem to be consistent with predictions from CGC initial 
conditions (IP-Glasma) without a hydrodynamic evolution, which may be key to understand the role 
between the initial condition and the collective expansion in such small systems. This is however 
still controversial, inasmuch as some hydrodynamic calculations are shown to reproduce the observed 
collective flow nature in p-Pb [315H317] . 


6 More about on Relativistic Hydrodynamics 

As discussed until now, the success of hydrodynamic modeling in nucleus-nucleus collisions is over¬ 
whelming and hydrodynamics is now one of the fundamental tools for the study of the initial state and 
transport properties of QCD matter in terms of experimentally observed collective flow parameters. In 
spite of these successes, one should be aware of the fact that there still exist several questions to be 
clarified within the hydrodynamic approach used in relativistic heavy ion collisions. 

First, the hydrodynamic fit to experimental data is not unique. There are still model-dependent 
systematic uncertainties. In a recent study, S. Pratt et al. [ 318] proposed a method to constrain in a 

“Reprinted with permission from K. Aamodt et al., Two-pion Bose-Einstein correlations in central Pb-Pb collisions 
at v /sVv = 2.76 TeV, Phys. Lett. B696 (2011) 328. Copyright © 2010 CERN. 
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global manner the EoS in hydrodynamic modeling using the conditions imposed by observables, such as 
particle spectra, HBT radii, and elliptic flow coefficients. In Fig. 11 we see that the overall shape of the 
EoS is considerably constrained but there is still a large class of different EoS can be equivalently good 
with respect to the observables selected. Of course, not all the constraints available from the data have 
been considered in this example, such as the higher order Fourier flow components, and the inclusion 
of additional information may constrain the EoS even further, which makes such kind of analysis very 
promising based on the state-of-art statistical methods. Bnt at the same time, in the present analysis, 



T (MeV) 

Figure 11: (color online) Equation of state calculations from lattice QCD and from models unconstrained 
(left panel) and constrained (right panel) by data from RHIC and LHC. Figur^] taken from [318]. 

no effects of viscosities, especially their temperature dependence, are taken into account. Also, several 
ambiguities which come from the initial condition and the final state hadronic cascade part are absent. 
Therefore, we still cannot conclude that the hydrodynamic representation is uniquely determined at the 
present stage. 

Second, according to the Little Bang scenario, in the central rapidity region a hot strongly coupled 
QGP is formed at a very early time after the collision (~10 -1 fm) in almost thermal equilibrium. 
Unfortunately we do not know yet what is the exact mechanism behind the thermalization process 
though several approaches have been proposed [3191 1324] , Furthermore, the success of the hydrodynamic 
description does not necessarily imply that the matter in question is in LTE. For a peculiar example, we 
may consider an ideal gas of massless particles. When the system is locally isotropic in the momentum 
space, then independently of the momentum distribution, the energy-momentum tensor becomes exactly 
that of an ideal fluid with the “EoS”, e = 3 p in LRF. See also Ref. [325] . 

In addition, at the present stage the great success of the hydrodynamic approach is mainly concen¬ 
trated in the mid-rapidity region of ultra-relativistic collisions where the baryonic chemical potential 
is close to zero, i.e., where the matter created from the vacuum (converted from the incident energy) 
is dominant. When we extend the hydrodynamic description to the non-central rapidity domain or to 
lower energies to describe the BES at RHIC, CBM at FAIR, and NICA programs, one needs a hydrody¬ 
namic model which accounts for the high baryonic current and the corresponding EoS at finite baryonic 
chemical potential. Furthermore, in the domain where the baryonic charge current becomes dominant 
we expect significant effects of non-equilibrium since these currents strongly preserve the information 
about the initial incident state. In order to discuss quantitatively these situations, one still has to clarify 
the validity and the resolution power of the hydrodynamics. 

“Reprinted with permission from S. Pratt et al, Constraining the Equation of State of Superhadronic Matter from 
Heavy-Ion Collisions , Phys. Rev. Lett. 114 (2015) 202301. Copyright © 2015 American Physical Society. 
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In this section we first discuss several approaches to relativistic viscous hydrodynamics and the 
meaning of the physical assumptions behind them. In particular, we discuss the concept of C-G hidden 
in the final form, which is essential to understand the resolution power of hydrodynamic approaches 
when we use them as a tool to extract the initial condition of energy and momentum distributions from 
the observed flow profiles. 


6.1 Different formulations of relativistic dissipative hydrodynamics 

There are several approaches used in the formulation of relativistic dissipative hydrodynamics and, 
differently from the ideal case, the resulting equations are not unique. Thus, the theory of dissipative 
hydrodynamics is not yet completely established as we will see below. 


6.1.1 Macroscopic approaches 


In section 3.5, we briefly discussed the phenomenological derivation employed by Israel and Stewart. In 
that case, the usual thermodynamics is extended to include II, 7r / “', and as additional thermodynamic 
variables. Once the number of thermodynamic variables is changed, various thermodynamic laws are 
affected. However, there is no established theory for such a generalization of thermodynamics to include 
irreversible currents. 

One example of such a theory is called extended irreversible thermodynamics (EIT) In this 

theory, the first law of thermodynamics near equilibrium is generalized to be 


TdS = dE + PdV - fidN 


^Ud{UV) - —n^din^V) - —u^d^V). 

Tn T n T k 


(47) 


By calculating the entropy production with this new equation as is done in the linear irreversible 
thermodynamics, one obtains evolution equations for viscosities. The results are the same as those in 


Eq. (33). The internal-variable theory (IVT) proposes another generalization of thermodynamics, but 
the essential structure of the derivation is the same as in EIT [327] , 

In the examples above, n, and are considered to be additional thermodynamic variables. 
However, this is not the only way to extend thermodynamics. As an example, let us consider the 
divergent-type theory (DTT) (32814330] where thermodynamics is obtained from 


TdS» = -iidN^ - u u d,T lw - T^sdA^ 5 , 


(48) 


where and A fluS are additional thermodynamic variables. By using these variables, the viscous tensor 
is expressed as, 

** = ir - A*»g; _ i . (49) 

The evolution equation of ^ is induced by the velocity gradient as is the case of the shear viscosity. 
Another phenomenological method which can be used to derive hydrodynamics is known as GENERIC. 
We refer the reader to the Refs. (33111332| for further details. 

There is an important property in phenomenological derivations which is known as Curie’s principle. 
In linear irreversible thermodynamics [98], irreversible currents with different tensor properties cannot 
be coupled. In hydrodynamics, we consider three thermodynamic forces, d^u u , and d^(n/T) but 
these terms are not consider to mix. For example, we can consider a vector current formed by a linear 
combination of thermodynamic forces such as 

F»d„u I ' + Gd^, (50) 
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where F M and G are the vector and scalar Onsager coefficients, respectively. This is obviously a 
mathematically covariant vector but its inclusion is not allowed from Curie’s principle, where any 
thermodynamic force of a given rank tensor can induce only irreversible currents of lower or equal 
rank. In non-relativistic fluids, this principle is experimentally established within linear irreversible 
thermodynamics. 

The applicability of Curie’s principle in relativistic and/or non-linear regimes is not clear, but 
assuming that it can be employed, one can reduce the number of terms which can possibly appear just 


based on Lorentz invariance. For example, we assume that Q M is given by Eq. (32) in the derivation of 


the Israel-Stewart theory [100] . Therefore, it is possible to extend the form of Q M to include the cross 
terms II and z/h In fact, such terms are considered in the original argument by Israel and Stewart. 
However, those terms lead to a coupling between n and which is then forbidden by Curie’s principle. 


6.1.2 Mesoscopic approaches (Truncation of Boltzmann equation) 

The relativistic Boltzmann equation for a rarefied non-degenerate gas of classical particle^ of mass m, 
without external forces, can be written as 


P 




(51) 


where the right hand side is called the collision term, which describes the changes of the one particle 
distribution function / (x,p) in terms of binary collisions, 


P 


- p ^ —> p ^ 


(52) 


with the differential cross section denoted by dcr/dVl. Note that f^' 1 = f and $ is the invariant 

flux, = y jF'p/j 1 —m 2 . In the collision term, energy and momentum are conserved and these con¬ 


servation laws reduce the integration of the hnal states specified by from six to two, which 

corresponds to the integral of the solid angle Q between p and pP\ Since the integrand of the collision 
term is Lorentz invariant, we can choose the center of mass system of a given p,pP> so that the final 
momenta are trivially expressed as functions of the solid angle. We also have to integrate over all the 
colliding momenta, p^ on p. 

The Boltzmann equation is an effective theory to describe the non-equilibrium dynamics of a rarefied 
gas. This is obtained by, at least in classical non-relativistic case, the lowest order truncation of 
the exact many-body dynamics called Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy 
equation [33511336] . That is, the two particle distribution function which appears in the collision term 
of the BBGKY equation is approximately replaced by the products of the one particle distribution 
function as is shown on the right hand side of Eq. (51). This truncation is called Stosszahlansatz or 


molecular chaos assumption and, due to this approximation, the time reversal symmetry is violated in 
the Boltzmann equation. This is a kind of coarse-graining procedure, leading to the existence of the so- 
called H-function, which decreases in time. Furthermore, such a truncation is consistent if the system 
is dilute. Therefore, the Boltzmann approach is expected to provide a solid theoretical framework 
to describe the transport properties involved in the dynamics of a relativistic dilute gas such as, for 
example, the hadron gas phase present in the kinetic freezeout stage of relativistic heavy ion collisions. 

On the other hand, since the microscopic interactions are considered as point-like (which is consistent 
with molecular chaos), any hysteresis effects in the dynamics are not taken into account. In fact, it 
is known that in the non-relativistic case the dynamic evolution of fluids described by the Boltzmann 
equation and by the Navier-Stokes-Fourier (NSF) theory show some qualitative differences. The famous 


3 Effects of quantum statistics according to fermion or boson can also be treated 1333) . 
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example is the velocity correlation functions. The value obtained from the non-relativistic Boltzmann 
equation decays exponentially in time [1337 ] while that from NSF theory one finds a power law tail 
[ 335; . 336 ;. 13381 - 1340] . In addition, the bulk viscosity identically vanishes in the Boltzmann approach 
for the non-relativistic case. Therefore, one has to be aware of the applicability and limitations of the 
Boltzmann approach. 

Nevertheless, when extracting the long-time-scale dynamics described by the linearized Boltzmann 
equation it is known that one can reproduce the form of NSF theory. That is, the linearized Boltzmann 
equation transfers some basic features of the microscopic collisional processes to the macroscopic viscous 
dynamics. Because of this, it is important to investigate how hydrodynamics emerges from the Boltz¬ 
mann equation especially in the relativistic domain where one does not have a firm phenomenological 
guidance as shown in the previous subsection. 

The Chapman-Enskog method [ 334 . 336 ] is a typical method employed to obtain macroscopic equa¬ 
tions of motion from the Boltzmann equation. As mentioned above, the Boltzmann equation itself is 
already an effective theory, but it still contains rapid motions that are irrelevant in the hydrodynamic 
regime. Therefore, in this case we are not interested in the exact behavior of the solution and, in fact, 
we want to extract an effective one-particle distribution function of the type f(x,p ) = fo(x,p ) + 5f 
which characterizes the long-time-scale dynamics of the relativistic Boltzmann equation, where / 0 and 
5f are defined below. 

How can one determine the macroscopic dynamics? In quantum mechanics, the time evolution of 
the wave function is expressed in terms of a linear combination of the energy eigenfunctions. Then, the 
scale of time oscillations is characterized by the magnitude of the energy eigenvalues. Here we want 
to apply a similar argument to proceed with the coarse-graining reduction scheme of the Boltzmann 
equation. For this, we consider the following linearized Boltzmann equation, 

p>*dJ = foC<f>, (53) 


where 0 = 5f/fo = f /fo — 1. Here we introduced the linear collision operator C defined by 

c*= [ /« [(tfW + - (<t> m + 


do 4 :i //' J 


(54) 


Po 


In the derivation of the collision operator above, we required that the basis function /o should satisfy 
the detailed balance relation (see below), 

fo (. V ) fo (P (1) ) = fo (p (2) ) fo ( P (3) ) , (55) 

or equivalently 

ln/o + ln/ 0 (1) =ln/® + ln/ 0 (3) . (56) 

The eigenvalue problem of this operator plays a crucial role. The collision term of the linearized 
Boltzmann equation contains various processes whose scales are characterized by the magnitude of the 
eigenvalues of the collision operator, 

C (ft = A0“, (57) 

where a is the degeneracy index of the eigenvalue A. Any function 0 which describes the changes of 
/ due to the effect of binary collisions should be expressed in terms of a linear combination of these 
eigenfunctions, (f> = ^ A a c A , a 0“. 

It is known that all the eigenvalues are real and non-positive [334] . The slowest modes have zero 


eigenvalue A = 0. It is obvious from Eq. (54) that any function 0 (p) which is conserved in the collision 
process Eq. ( f52| is an eigenfunction of C with a null eigenvalue. In a binary collision, as is well- 
known, there are five conserved quantities: particle number, energy, and momentum. Thus, they are 
degenerated eigenfunctions with A = 0, 


Cl = C p» = 0. 


(58) 
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On the other hand, due to the condition Eq. (56), In f 0 is an eigenfunction of C with A = 0. Since 


there are no other independent conserved quantities than Eq. (58), the most general form of f 0 should 
be expressed as 

In f 0 = a(x)l- ^{x)p^. (59) 

As was pointed out before, when there are conserved quantities the time scales of the corresponding 
densities are very long compared to the microscopic time scale. Therefore, it is natural to choose such 
/o as the starting point of the C-G reduction procedure with respect to their time-scales. Additionally, 
the above /o is nothing but the Maxwell-Jiittner distribution and this expansion is equivalent to the 
expansion around the equilibrium distribution function. 

Now, suppose that the amplitude of </> is very small (\(j)\ <C 1) so that the effective distribution is 
essentially given by / = /o- Then one can show that the dynamics described by Eq. (53) is reduced to 
ideal relativistic hydrodynamics by identifying 


T^(x) = (p^) Uo) , 

where we defined the notation for the average of any function O ( p ) by /o as 

1 


(°Xfo) ~ 


(2tt) 


(p) f 0 (x,p). 


Using the Landau definition, Eq. (22) for the fluid velocity tG, one can in fact show that 


T» v (x) = (e + P) uW - g^P, 


with 


£ = <Kp") 2 )(a). 


p = --{a '“'j W ,) (A) . 


(60) 


(61) 


(62) 


(63) 


Note that, for consistency, the time scale of the variation in a(x) and /3 M (x) should be much slower than 
the microscopic ones. Because In f 0 describes the motion in the functional space spanned only by the 
five eigenfunctions with A = 0 the dynamics of the parameters a and are completely disentangled 
from the change in the momentum states due to the collision term. Therefore, within the time scale 
corresponding to the change of a and (3^, detailed balance is considered to be held. That is, the evolution 
only through the changes in a and (3^ corresponds to those of the quasi-static process that maintains 
the local equilibrium state of the gas. In other words, the dynamics is time-reversible, i.e., that of an 
ideal fluid. 

Therefore, to consider the violation of the detailed balance, which is directly related to the effect of 
dissipation, we need to expand the dynamical functional space by introducing other eigenfunctions of 
C with finite eigenvalues. In the Chapman-Enskog method, this extension is implemented by assuming 
that the Erst order deviation from the space of A = 0 is described by the linearized Boltzmann equation 
(53) (see Ref. [332]). By definition, the deviation of the one-particle distribution should relax at scales 


quicker than the time scale of /o and thus, in the long-time limit, the right hand side of Eq. (53) is 
approximately replaced by 

p^dj « p^djo- (64) 

Substituting this into the linearized Boltzmann equation, the deviation is given formally by 


- c-ya, ; in/„. 


(65) 


Now the right-hand side contains a component proportional to p^p v which is not a collision invariant 
anymore and the deviation from /o is affected by a component which involves a faster time scale 
than those from the collision invariants. Substituting this into the linearized Boltzmann equation, we 
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obtain the relativistic NSF theory. This is the physics behind the coarse-graining of time scales in 
the Chapman-Enskog method and the above-mentioned scheme is mathematically implemented as a 
systematic expansion in powers of the Knudsen number. 

Relevant contributions which survive even after taking the long-time limit are sometimes called 
secular terms in the derivation of coarse-grained dynamics. For example, it is known that the systematic 
collection of the secular terms is important to obtain a master equation in other mechanical systems 
such as celestial mechanics and quantum mechanics [ 341. 1342] . The C-G procedure employed here 
corresponds to the van Hove limit. However, the way the secular terms can be resummed is neither 
trivial nor unique. In Ref. [ 343j , for example, the method to reduce nonlinear differential equations is 
applied to the relativistic Boltzmann equation and another relativistic generalization of the NSF theory 
is obtained. 

In the Chapman-Enskog method the perturbative contribution to the effective distribution is ex¬ 
pressed in terms of / 0 which is a function of the collisional invariants. On the other hand, to obtain a 
causal model, one needs to introduce new hydrodynamic variables as mentioned before. Naively think¬ 
ing, such additional terms can be obtained by just generalizing the form of /o itself. However, this is 
a very difficult task and this constitutes one of the reasons why it is not trivial at all to obtain causal 
hydrodynamic models using the Chapman-Enskog approach. 

The moments method, which was proposed by Grad [ 344'j . is a representative method to obtain 
causal hydrodynamic models from the Boltzmann equation. In this approach, the Boltzmann equation 
is converted to coupled equations for the moments, which are, for example, defined by 

MT 2 '= f = (Kj’Ti> < “‘-'P“ ,> }w), (66) 

where p <otl ■ ■ ■ p ai> = P' 1 ■ • • p' 1 with being the symmetric traceless projection operator 

for any of two cq’s or two 00s [345] . One can obtain the exact coupled equations for these moments, 
which is equivalent to the Boltzmann equation if the expansion converges. 

In the Chapman-Enskog method, 0 (or equivalently 5f) is determined dynamically by solving the 
linearized Boltzmann equation. On the other hand, in the moments method implemented by Israel and 
Stewart, a functional ansatz of 0 is introduced, 


0 = Ain + 


(67) 


where 


n = --(A V" = (A '■"“VwW V = <?%/>■ 


( 68 ) 


This ansatz for 0 is called 14 moment approximation [ 100] , By assumption, the time evolution of 0 is 
restricted in the dynamical functional space spanned by n, z/ M , and 7r MI/ . As in the case of Chapman- 
Enskog’s method, one assumes that the deviation from /o is described by the linearized Boltzmann 


equation (53). Note that the above assumption for 0 resembles the Kawasaki relation which gives the 
relation between equilibrium and non-equilibrium distribution functions [346]. Once dynamical equa¬ 


tions for these new variables are obtained in a closed form, one can derive the dissipative hydrodynamics 
equations. 

However, closing the system of dynamical equations is not granted. In the moments method, one 
introduces a new ambiguity in the choice of the moments to derive dissipative hydrodynamics which does 
not have a counterpart in the Chapman-Enskog method [337] . For example, in the original argument 
of Israel and Stewart, the equation for the bulk viscosity is obtained from the moment equation of 
UvU\d fl (p fl p u p x ) It is however possible to obtain the equation for the bulk viscosity from another 
moment equation, for example, A4o- The form of the transport coefficients calculated in the moments 
method depends on the choice of the equation for the moments, which are used to close the coarse¬ 
grained dynamics [348] . 
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It is possible to close the moment equation without using the 14 moment approximation (67). In 


this approach, the form of the effective particle distribution function / is not specifically calculated. 
Instead, the linearized Boltzmann equation is mapped into a set of infinitely coupled equations for 
the moments and a coarse-graining reduction with respect to the time scale is applied directly to the 
moment equations by introducing an expansion around a diagonalized basis of moments of the collision 
operator [333] . 

The equations for the moments A4 r , A4 S , ■ ■ ■ have contributions from the collision operator. The 
corresponding terms can be expanded as 


JV, 


C«i w = / dT r (u^y pP <>“ ■ ■ ■r ,> foC<l> = - 4 2xi 


Ml —Mi 


(69) 


71=0 


The above equation is equivalent to the Boltzmann equation. To obtain C-G dynamics, as is done in the 
Chapman-Enskog expansion, where the eigenfunctions of the collision operator with the null eigenvalue 
A = 0 play an important role, one needs to introduce an appropriate truncation scheme in the system of 
moments. For the derivation of the hydrodynamic equation, the moments for l — 0,1,2 are important. 
Each Ni is in principle oo, but we introduce the truncation by choosing hnite numbers No, N\ and N 2 
as three parameters of this approach. Note that AR and AR for l = 0 and for l — 1 are excluded 
from the basis of the expansion. 

I 11 the basis where the matrix Ar\ is diagonal, the moment AR r " w is changed into A^ r " w . Let 
us assume that the deviation from equilibrium can be spanned only by A 0 , Aq and X^ u> among all 
possible moments Xft 1 "'** 1 . This corresponds to an expansion of the dynamical functional space in terms 
of N 0 — 2 scalar moments (A4 0 , AR ■ ■ ■ M-n 0 ), Ni — 1 vector moments (AR, AR , ■ ■ ■ A4R) and N 2 tensor 
moments (Mq V , ■ ■ ■ M^ 2 )■ 

The moments equations for A 0 , AR and Xq^ u> still couple to other higher order moments. Then 
the truncation of the coupled moment equation for AR ”' w is carried out by maintaining the first order 
terms of the expansion in terms of two parameters, the Knudsen number and the inverse Reynolds 
number. 

To express hydrodynamics in terms of II, zA and we change from the diagonalized basis to 
the original basis of AR r ” w . Then, for example, AR, which is associated with the bulk viscosity, is 
expressed as a linear combination of A 0 , A 3 ■ ■ ■ Xn 0 . To eliminate the dynamics of A "3 • • • Ajv 0 , we further 
assume that the time evolution of these moments is much faster than Ao and the values are replaced 
by the stationary solution of the corresponding equations for the moments. 

The hydrodynamic theory which is obtained in this scheme is called transient hydrodynamics. One 
can easily see that this hydrodynamic theory leads to a different model from the one obtained from 
the 14 moment approximation in general because the dynamical functional space used to perform the 
coarse-grained dynamics is different. In the 14 moment approximation, the dynamics of 0 is considered 
to take place in the space spanned by II, and 7 AO On the other hand, in transient hydrodynamics, 
0 is not written explicitly and, thus, a direct comparison is not trivial to be performed. However, 
the moment equations are expanded in the basis spanned by Ao — 2 scalar moments, N± — 1 vector 
moments, and IV 2 tensor moments and the scale of coarse-graining is different than that obtained in the 
14 moment approximation. 

Interestingly enough, if we choose Ao = 2, N\ = 1, and IV 2 = 0, then the moment equations 
are expanded in terms of Nio, A4q and in transient hydrodynamics, which is similar to the 14 
moment approximation. In fact, transient hydrodynamics leads to the same result as the 14 moment 
approximation in this case. There is also another approach based on the Boltzmann equation, cf. Ref. 

[335] . 

It is often claimed that hydrodynamics is obtained from an underlying microscopic theory by expand¬ 
ing in powers of the Knudsen number. This view seems to be based on the Chapman-Enskog method for 
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the Boltzmann equation discussed above. However, such a perspective is quite suspicious because one 
cannot improve upon the coarse-graining procedure even if one considers higher order correlations in 
the Chapman-Enskog method. In fact, non-relativistic hydrodynamics including such correction terms 
is called the Burnett and super Burnett equations [ 347] . which are known to suffer from an instability 
called Bobylev instability. Moreover, there are no examples of fluids which are correctly described by 
either the Burnett or super Burnett equations. Therefore, it is not obvious whether one should consider 
that C-G can be systematically implemented via the Knudsen number expansion. 

Finally, we have one more comment on the so-called derivation of hydrodynamics from the Boltz¬ 
mann equation. As was mentioned above, the Boltzmann equation is appropriate to describe dilute 
systems while the hydrodynamic description works better in dense systems. Therefore, the discussion 
presented here should not be interpreted as a rigorous derivation of the hydrodynamic equations from 
the Boltzmann equation itself. Rather, it is expected that the asymptotic behavior towards equilibrium 
in the Boltzmann equation correctly describes the deviation from a given LTE and in fact it can re¬ 
produce the structure of the NSF theory. This is because the asymptotic LTE state is an input in the 
case of linearized Boltzmann equation. Furthermore, we note that the dynamics described by the exact 
Boltzmann equation does not coincide with NSF theory as mentioned before. 


6.1.3 Projection Operator Method 

Microscopic dynamics has time-reversal symmetry and special care is need to violate it to obtain C-G 
dissipative dynamics. In this subsection, we explain how this is done in the framework of the projection 
operator method. 

The most traditional way to implement C-G in microscopic dynamics is the projection operator 
method in statistical physics [102 , 133511350H352] . In this approach, we first choose the dynamical space 
onto which our non-equilibrium dynamics is mapped and define a projection operator P to implement 
the coarse-graining with the chosen variables. By using this projection, we can systematically project 
out irrelevant rapid motions contained in a microscopic dynamics. 

Let us define the (time-independent) projection operator satisfying the following general properties, 


P 2 = P, Q = l-p. (70) 

From these definitions, one can see that Q extracts the dynamics orthogonal to P, that is, the irrelevant 
rapidly moving degrees of freedom. By using the properties of the projection operators mentioned above, 
the Heisenberg equation of motion for an arbitrary operator O is expressed as 

d t O(t) = e iLt PiLO + [ dse iL ^~ s) PiLQe iLQs iLO + Qe iLQt iLO , (71) 

Jo 

where L denotes the Liouvillc operator defined by 


iLO 


i[H,0\ for quantum systems, 
{H, 0} PB for classical systems. 


(72) 


Here {H, 0} PB is the Poisson bracket of H and O. This equation is called time-convolution equation 
of the projection operator and can be derived without specifying the concrete definition of P. The first 
term on the right-hand side is called the free-streaming term and corresponds to a collective motion 
which does not produce entropy. The second term represents the dissipative part of the time evolution 
while the third term is interpreted as the noise term. In the following, we exclusively discuss quantum 
mechanical systems. 

There are several possible choices for P, but we consider here the so-called Mori’s projection operator 
in the following discussion. Let us consider a set of macroscopic variables by an n-dimensional vector, 
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A 7 = {A]_, ■ ■ ■ , A n ). Then the Mori projection operator is defined as 


PO = Y J c i A il (73) 

i— 1 

where the coefficient is 

n 

Ci = J2(0,A])(A,A^. (74) 

3 = 1 

The inner product is given by Kubo’s canonical correlation, which is defined by 

(X, Y) = J* jTi [p eq e XH Xe- m Y] = j* j(e XH Xe- m Y) eq , (75) 

where p eq = e _/ 3 J 7 /Tr[e _/37/ ] with j3 being the inverse temperature. See Refs. 13531 13571 for details. 

Let us consider some type of microscopic dynamics characterized by the Hamiltonian H and derive 
its coarse-grained dynamics assuming that the dynamics is described by only one macroscopic variable, 
say, B. Then the macroscopic evolution equation of this macroscopic variable B is obtained from Eq. 
( ffTj ) as 

| ] - t B (*) = dsD(s)B(t - s) + f (t), (76) 

where 

= £(t) — Qt iLQt iLB. (77) 

Here we assumed that B is an Hermitian operator. This equation is still equivalent to the Heisenberg 
equation of motion and, to obtain dissipative behavior, we need to introduce the coarse-graining of the 
time scale. 

Note that if the projection operator P is chosen appropriately in order to completely extract macro¬ 
scopic motions, the time scale of the noise term £(t) is very short because it is proportional to Q and 
we can ignore the characteristic time scale of (£(£),£(0)) in D(t). Then we can assume that the time 
dependence of D(t) is approximately proportional to 5(t) and the exact equation above is reduced to 




(78) 


where 7 = J Q ^° dsD(s). This is called the time-convolutionless limit (TCL), which is equivalent to a 
Markov approximation [355 . 356] . We shall soon see a concrete example of this equation. As seen in 
Eq. (78) the equation obtained in the projection operator method contains a noise term that is related 
to the viscosity through the fluctuation-dissipation theorem. The relativistic hydrodynamic model with 
noise terms is discussed in Ref. [358] and references therein. However, this type of approach is known 
to present difficulties such as Lorentz covariance and stability [359] . 

The choice of the hydrodynamic variables is closely related to the applicability of the TCL, that 
is, the consistent violation of the time reversal symmetry. One can easily see that the TCL is not 
applicable when the time correlation function of the noise converges to a finite value in the long time 
limit, 

lim D{t) = const , (79) 

t—> OO 


because then 7 diverges. This means that macroscopic dynamics has not yet been extracted and we 
need to extend the definition of the projection operator P. 
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For example, in dynamics of diffusion, it is known that the projection operator is defined with two 
macroscopic variables, n(x) and J(x), which satisfy the continuity equation, to apply the TCL limitQ 
d t n(x) = —V ■ J(x). Then the equation corresponding to Eq. (78) is given by 


J^i(x,t) + VJ(x, t) = 0, 
d 

t d — J(x,t) + J(x,t) = -DVn(x,f), 

where the diffusion coefficient D and the relaxation time To are, respectively, given by 


dn 


D = — ( — ) lim lim Gj (cu, k), — = - 

T 


D lEJA(f(x),J-(0)) 


dnjrr^-t ok^o J ' ’ td 3 f d 3 x(3n(x), 3n( 0)) 


where SA = A — {A) eq and the retarded Green function is 


(80) 

(81) 

(82) 


Gj(u,k) = / dt 


^ xe —V :kx ([J ?: (x,t),J'(0,0)]) eq . 


(83) 


Here the noise term has been dropped. 

The projection operator approach provides information on how one should violate the time reversal 
symmetry. For example, one may consider that the TCL is the lowest order truncation of the following 
expression of the time-convolution integral term, 


dsD(s)B(t 


s = 


OO 

£ 

72—0 ' 


dsD(s 


d n B(t ) 

'“a v 


(84) 


and consider that the solution can be improved by considering higher order terms but that is not 
true. First of all, the TCL is an artificial operation to violate the time reversal symmetry and this 
is done not for the purpose of the simplification of dynamics. In fact, it is known that the TCL is 
the most appropriate procedure for the violation of time-reversal symmetry in the projection operator 
method because it can pick up all the secular terms correctly and an exact relation, called f-sum rule, 
is maintained [356]. Alternatively, if we consider a correction term to the TCL, secular and non-secular 
terms are considered on an equal footing and the f-sum rule is violated [ 356 , 35? j. If we want to take into 
account more microscopic dynamics associated with D(t), we need to extend the number of dynamical 
variables to define the projection operator as was done in the example above. 

The problem of the projection operator method is that it is difficult to maintain the manifest 
Lorentz covariance. The attempt to formulate relativistic hydrodynamics is, thus, so far limited to 
the derivation of linearized equation [361] . It is however still possible to utilize it to define transport 
coefficients because these are Lorentz scalar quantities and, in this case, it is enough to consider linear 
deviations from equilibrium. In fact, it is known that this approach leads to the correct transport 
coefficients that appear in NSF theory [ 354 j. 

To obtain, for example, the shear viscosity in a causal hydrodynamic model the projection space is 
defined by the two variables, A = T yx and T 0x . Then r) and the corresponding relaxation time r n are 
defined by 

V = _ Vgkn _ V = f d 3 x(T yx (x),T yx ( 0)) 

/ 3{e + P ) P 2 J d 3 x(T 0x (x) 1 T 0x (0)) , T n (e + P) f d 3 x(T 0x (x),T 0x (0))' [ J 

4 In the textbooks of statistical physics [335(1350 . 360j, it is written that the diffusion equation can be obtained in the 
projection operator method. However, in the derivation, an artificial approximation is introduced. If this approximation 
is not used, we encounter the singular behavior of the diffusion coefficient and the definition of the projection operator 
should be generalized. This fact was first pointed out in Ref. 13571 . 
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where i~Igkn is the usual result of the shear viscosity coefficient in the relativistic NSF theory [102113551 
[362] , This expression is called Green-Kubo-Nakano (or simply Kubo) formula, but, exactly speaking, 
this result was obtained using Zubarev’s non-equilibrium statistical operator method [363] . Note that 
the results above are universal and are thus valid not only at finite temperature but also at finite 
temperature and nonzero density |356| . The results for the bulk viscosity coefficient are shown in Refs. 
[1021 13551 1356] . 


The shear viscosity result has an interesting relation to other kinetic theory calculations and also 
AdS/CFT calculations. In fact, the lowest order calculation of Eq. (85) without including the pair 
creation-annihilation effect gives the exact same result as the moments method within the 14 moment 
approximation obtained in Ref. 136 1 . However, the contributions from the pair creation-annihilation 
effect can be as large as 20 % in Eq. (85) and 80 % for the corresponding quantity for the bulk viscosity. 
Thus, one cannot ignore this quantum effect. 


The transport coefficients of transient hydrodynamics assume different values as the number of 
moments in the expansion is increased. However, this difference is reasonable because the expansion 
basis of the dynamical functional space is now different. In the definition of the projection operator 
above, the deviation from thermal equilibrium is expressed by the T 0x and T yx , the former corresponds 
to (e + P)u x and the latter 7i yx . This corresponds to expand the dynamical functional space using the 14 
moments as was done in the moments method for the Boltzmann equation. To reproduce the result of, 
for example, the 23 moments approximation done in transient hydrodynamics, one needs to expand the 
definition of the projection operator. For such a systematic generalization of the projection operator, 
see Ref. [353] . 


On the other hand, if we consider the next order correction of the expansion by Eq. (84), the form 
of Eq. (85) is modified as was pointed out in Eq. (65) of Ref. 1 355 . 356] . and it coincides with the same 
expression obtained from AdS/CFT. See also Ref. [(352] where it was shown how one can reproduce the 


result of 14 moments approximation by modifying the AdS/CFT calculation. 


It is also worth mentioning that we need to use the modified definition of the energy-momentum 
tensor for the scalar held theory in calculating the bulk viscosity consistent with conformal symmetry 
[3551. 356]. It should also be noted that there are two philosophically different concepts involving the 
transport coefficients. When an irreversible current is induced by the application of an external held, as 
is the case of electric conductivity, such a perturbation can be expressed in the form of the interaction 
Hamiltonian. On the other hand, phenomena like heat conduction and viscosity are induced by the 
changes of the boundary conditions and, thus, these effects cannot be expressed as the interaction term. 
The importance of this difference was already pointed out by Kubo himself. See also the discussions in 
Ref. [366] . 

As we have discussed, there are two main factors needed to obtain the C-G dynamics from a 
given microscopic theory; one is the choice of macroscopic variables and the other is the Markov limit. 
These are, however, not independent. To employ the Markov limit, we need to choose a complete set of 
macroscopic variables. The difficulty of the C-G procedure is attributed to the absence of the established 
criterion to choose macroscopic variables, which is a common problem not only in the macroscopic and 
mesoscopic derivations which we have discussed, but also in the derivation of conformal fluids which will 
be discussed soon later. In the present case we chose conserved densities and the corresponding currents 
as macroscopic variables, but there are situations where one may have to extend this definition. For 
example, if there is a continuous phase transition, because of the critical slowing down the time scale of 
the order parameter held becomes sufficiently slow and we need to consider it as another macroscopic 
variable [367] . 
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6.1.4 Conformal hydrodynamics and transport coefficients with AdS/CFT correspon¬ 
dence 

Another method to derive relativistic hydrodynamics from the microscopic degrees of freedom in the 
strong coupling limit has been developed within the framework of AdS/CFT correspondence. There are 
various studies involving conformal hydrodynamics |368j (in which the hydrodynamic equations change 
covariantly under Weyl transforms of the metric) and recently this approach has been extended to study 
nonconformal theories as well [369] . The derivation of hydrodynamics itself is phenomenological, but one 
can calculate the transport coefficients from a microscopic theory using the AdS/CFT correspondence, 
which means the equivalence between a conformal held theory such as Af = 4 supersymmetric Yang- 
Mills theory and superstring theory in anti-de Sitter spaces [368] . 

First we discuss the derivation of hydrodynamics in this case. The philosophy of this approach 
is based on the interpretation that hydrodynamics can be obtained from kinetic theory within the 
Chapman-Enskog method which can be regarded as the Knudsen number (or equivalently derivative) 
expansion of the linearized Boltzmann equation as we have discussed. As it has been discussed so often 
in other approaches, one first specifies the hydrodynamic variables. We choose slow variables associated 
with the conservation of energy and momentum, e and tffi. If hydrodynamics can be expressed only 
by those variables, one obtains the ideal hydrodynamic case which we have already discussed. To 
obtain a dissipative theory, one needs to include the deviation from equilibrium. In this approach, 
this is assumed to be given as a function of e and d fJj u 1 ' and consider up to the second order terms in 
derivatives. Another method used in this derivation can be found in Ref. [370]. 

The advantage of considering a conformal fluid is that the continuity equation of the energy- 
momentum tensor should transform covariantly under Weyl transformations and, thus, the number 
of second order terms are reduced. 

The hydrodynamic theory derived this way has a term associated with the vorticity, 

5T- = Ia'-A^cU/? - dpu a ). (86) 

Such a term does not appear in the kinetic derivation of the Chapman-Enskog method and the 14 
moment approximation of the moments method and, thus, has been considered as a quantum effect. 
However, it was found that the same term can be naturally obtained in kinetic theory using transient 
hydrodynamics |345j . 

One can easily see that this is similar to the second order correction in the Chapman-Enskog method 
and then what we obtain is the relativistic generalization of the Burnett equation. In fact, it is known 
that this equation contains an unphysical propagation mode [368]. To obtain a causal model, one needs 
to replace some of the terms in the equation by the time derivative of by hand. 

Another advantage of discussing a conformal fluid is that one can compute the transport coefficients 
directly from the AdS/CFT correspondence. The retarded Green’s function of the energy-momentum 
tensor can be calculated from two different ways; the Erst one involves the energy momentum tensor 
equation which was derived phenomenologically above while the other one involves the microscopic 
theory, such as supersymmetric Yang-Mills theory, via the AdS/CFT correspondence. Expressing the 
retarded Green’s function in the Fourier variables uj and k and comparing both results at the same 
order one can find the desired transport coefficients. This method is essentially the same as the indirect 
Kubo method Ref. [366]. For example, the shear viscosity coefficient obtained is given by the famous 
ratio 

r) 1 

S 47T 

This result was first obtained in Refs. [ 157. 371 ] and it was considered to be a lower bound for the value 
of the shear viscosity coefficient, called the KSS bound. 
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There are several attempts to generalize the approach mentioned above to non-conformal fluids, see 
Refs. 1372Vt373] . 

As was mentioned above, this derivation of hydrodynamics is inspired by the success of the kinetic 
derivation. However, as was also emphasized before, there is a critical difference between the kinetic 
derivation and the present microscopic derivation. That is, in the Boltzmann equation time reversal 
symmetry is already broken but the microscopic theory is still symmetric. In the above argument 
involving the AdS/CFT correspondence, it is implicitly assumed that the symmetry is appropriately 
violated by the truncation of the derivative expansion. The validity of the assumption should be 
carefully investigated. 

6.2 Variational Principle 

In this subsection we present a quite different view from what has been discussed until now. As 
mentioned before in this text, hydrodynamics is a classical effective theory for interacting matter, 
which emerges from the underlying microscopic degrees of freedom through some kind of coarse-graining 
procedure. There are many ways to introduce a given C-G scheme. In order to obtain a closed system of 
dynamical equations only in terms of these C-G hydrodynamic variables, one usually needs to introduce 
further some approximations or truncations. As discussed in the previous section, depending on how 
one introduces these truncations, many different forms of dissipative hydrodynamics can be obtained. 

A different way to obtain a closed system of dynamic equations for a set of given C-G variables is 
the variational method. It has been shown in the very early days that relativistic hydrodynamics can 
also be formulated in terms of a variational principle [ 375. 376] . Thus, if we express the model action in 
terms of the selected hydrodynamic variables, the variational principle leads to the optimal dynamics 
for these variables. Below we show how the effective dynamics of C-G variables can be described in 
variational form. 


6.2.1 Coarse Graining in Variational Principle 


In order to perform the C-G explicitly, let us consider for example a classical microscopic system which 
contains a large number of quickly moving point-like particles. We can define the density 

pM) = (*)) 

i 

where the sum is done over all the particles in the system and r* (t) refers to the position of 7-th particle 
at time t. The symbol ~ is used to show that the variable is spiky (delta-function like) in space and 
rapidly changing in time. The corresponding current is 

j M) = J2^jr s ( r " ri (*))> 

i 

and we can easily see that they satisfy the continuity equation, 

dtp + V ■ j = 0. (88) 


However, we usually do not require a very precise resolution both in space and in time to describe 
collective flow behavior. Thus, we introduce an averaged smooth density distribution p (r, t) and current 
j from the original distributions using a four-dimensional smoothing kernel |377] as 


j (r, t) = j it' J dh'U ((' - t) W (r - r') j (r', t'), 


(89) 

(90) 
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where 


(91) 


and 



dt U (t) 


d 3 r W (r) = 1, 


U(t),W( r) -> 0, |£|,|r|»r, h 


(92) 


with r and h are given constants which characterize the time and space C-G scale, respectively. We 
can take, for example, U and W as Gaussian distributions with width r and h. It is easy to verify that 
dp/dt + V ■ j = 0. 

In an analogous manner, we can construct the smoothed (coarse-grained) energy-momentum tensor 
T M1/ as a convolution of the original very spiky and quickly changing T ^ u . It is easy to show that such 
an energy-momentum tensor again satisfies the continuity equation, = 0, as far as the original 

T , J,J does. Through the above C-G procedure, the two macroscopic fields, j M (r,t), and T M ^(r,t) can 
be constructed and, by definition, they are smooth and mildly changing within the scales of W and 

U. Now, instead of constructing the equation of motion from the microscopic dynamics for the original 
variables (in the present example, the set of large number of trajectories {rj (t)} ), we ask what is the 
optimal time evolution for j M (r,t), and (r, i), given that we know the initial and final values of 
them. 

Usually, the variational principle describes time-reversal symmetric problems. Thus, we first consider 
the ideal fluid case. In this case, the number of independent variables are reduced to five and they are 
specified by the velocity field u^ with = 1 and the two proper scalar densities, £ and n. They 
are related to j' 1 (r, £), and (r ,t) as n — and tG and £ are determined by the definition of 

Landau frame as the eigenvector and eigenvalue of (r, £). The number density n* in an observing 

frame is given by 717 , where 7 = u° is the Lorentz factor. Now, one needs to determine the form of the 
functional to be optimized. 


6.2.2 Hydrodynamic Action and the Ideal Fluid 


Now we would like to construct the variational principle to describe the dynamics of these C-G variables. 
By construction, the (0, 0) and (0, i ), i — 1, 2, 3 components of energy-momentum tensor T^ u represent 
the energy and the momentum densities carried by the fluid element, respectively. To stress this, let us 
denote them as 


H 
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i 


(93) 


Since E = f d 3 r'H is the total energy of the system, we can identify TL with the Hamiltonian density. 
Now, from Landau’s definition of the local rest frame T 3u u v = eu° we obtain 

—£ = Ti ■ v — H, (94) 

where v is the velocity of the fluid element. This equation indicates that the quantity — e should be 
the Lagrangian density of the system (L = pq — H). Therefore, the optimal dynamics will be given by 
modeling the action, 


/ = 


d A x e ( x ). 


(95) 


in terms of hydrodynamic variables, i.e., the local proper thermodynamic variables and the velocity 
held v. For example, for an adiabatic motion of fluid element, we can specify e as a function of any of 
the conserved proper densities, say 43 In terms of the density in the observational system of reference 
p — p */7 and the action is 


I {P , v] = 


d^x £ ( — 

7 


(96) 


D For instance, p can be chosen to be the entropy density, s in ideal hydrodynamics 
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where 7 = l/\/l — v 2 is the usual Lorentz factor. The variation should be taken with respect to 
the density p(x) and the velocity held v (x). In fact, it is known that this procedure reproduces the 
equations of motion ideal relativistic hydrodynamics [378] . 

Usually the variational procedure is performed in Eulerian coordinates and a constraint among p 
and v has been introduced using the Lagrange multiplier method. Below we sketch how to directly 


derive the equation of motion from the action Eq. (96) using Lagrange coordinates. 


We consider the situation where there is no chaotic motion and the time extension is finite so that 
we can follow the fluid motion for each infinitesimal fluid elements from the beginning. Let tr (t) be 
the space coordinate of a fluid element at time t whose initial position was R. We assume that tr (t) 
and R have an one to one correspondence for any t. The fluid velocity is vr (t) = dr^ (t) /dt. Let 
Po (R-) be the initial distribution of the fluid elements. Without any loss of generality, we can always 
choose the distribution of initial fluid elements to be p 0 — const. Then, the fluid density distribution 
p (r, t) is determined as 

p*(r,*)=jrf (97) 


where J = det |<9r/<9R|, is the Jacobian of the transformation R —» r = tr (t). Equation (97) is a direct 
consequence of the conservation of the number of particles in the fluid element, p(r,t)c/ 3 r = po<J 3 R. 


The action (96) can be re-expressed as 


I [{r R , Vr}] = j dt J d 3 R C [J, 7 ], 


with the Lagrangian density C 


C[J,i\ = ~Je 


Po 
7 J 


(98) 


(99) 


Note that the Lagrangian above contains only t and R derivatives of the dynamical variables tr (t) 
through the Lorentz factor 7 and the Jacobian J. Using the properties of the determinant, the term 
which appears in the Euler-Lagrange equation can be calculated as 


E 

k =1 


d 
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dC 


dJ 


dJ d(dr R /dR k 


T __ :ac 
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Furthermore, 
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pl8(l/p) P,> \p) 



( 100 ) 


where V = 1/p is the specific volume, U = e/p, is the specific energy, and P is the pressure. On the 
other hand, we find that 


dC 

dw 


(p1 

p dp dv \ 7 


* £ + P 

= Po —V. 


( 101 ) 


Thus, the action given in Eq. (98) leads exactly to Eq. (21), 
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dr 


e + P 
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( 102 ) 


Here we have shown that the relativistic Euler equation can be obtained from the variational prin¬ 


ciple, where the form of the action (96) was inferred from the Landau definition of the fluid velocity as 


the time-like eigenvalue of the energy-momentum tensor. However, once given the action Eq. (96), the 
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energy-momentum tensor can be derived as the Noether conserved current associated with the local 
translational invariance of the Lagrangian density [375] . 

Finally, we would like to make a comment on the applications of the variational approach on col¬ 
lective flow observables in relativistic heavy ion collisions. In this formalism, we specify the final state 
of the collision in terms of macroscopic hydrodynamic variables which correspond to these observables, 
and we also introduce some model initial condition consistent with the centrality event selection. Ob¬ 
viously, for a given event, the C-G variables are uniquely defined but the inverse is not true. There 
are many events with distinct microscopic configurations that might give the same hydrodynamic final 
state. This is also the case for the initial condition. Furthermore, the C-G procedure can even be 
extended to include the averaging procedure over different events which are rather loosely classified 
(see the discussion in Ref. [364] ). This means that, in practice we may further allow a large statistical 
ensemble contained in the C-G scheme within a certain given resolution of the physical definition of the 
hydrodynamic event. If the ensemble of all the microscopic dynamics (including the different events) 
equivalent within the C-G criteria is large enough, then the macroscopic variables specified by such C-G 
procedure distribute sharply around their mean-values as a consequence of the central limit theorem. 
The variational approach with these C-G variables then results in the hydrodynamic picture of the 
event averaged flow. This can be used to explain why smooth initial conditions reproduce the collective 
behavior of the event averaged collective observables reasonably well as the EbyE fluctuating initial 
condition hydrodynamics, and also why collective signals are so robust. 

6.2.3 Application of Variational Approach - Effective Hydrodynamics 

Once the variation principle is established for relativistic hydrodynamics, we can apply such an approach 
to obtain physical insight about the collective behavior by introducing a simplified parametrization of 
the fluid dynamics and approximately obtain the dynamical behavior of the solution. For example, the 
equation for the dynamics of a cavitation bubble, known as Rayleigh-Plesst equation and its relativistic 
form can be derived from the variational approach [578] . Also, many years ago in low energy nuclear 
physics the description of nuclear collective motions was often formulated in terms of the variational 
principle. For instance, the dynamical mixture of Steinwedel-Jensen and Gamow-Teller modes in giant 
resonance phenomena is discussed using the variational approach and it described the behavior of the 
observed data very well [580]. Such models can be considered as examples of the C-G procedure to 
extract the global collective dynamics. 

SPH Formalism The same line of thinking can also be applied to obtain a numerical method to 
solve the hydrodynamic equations, called Smoothed Particle Method (SPH). The SPH algorithm was 
first introduced in astrophysical applications [3811 1382 ) and in Refs. [12], 133] this numerical method 
was applied in relativistic heavy-ion collisions using the variational approach discussed in the preceding 
subsection. The basic idea of the SPH method is to parametrize the continuous density distribution 
of any extensive physical quantity in terms of a discrete sum of base functions. This is somewhat the 
inverse way of the coarse-graining procedure for particle systems. Let a(r, t) be the density distribution 
of a conserved extensive quantity. Then we parametrize 

N 

a(r,t) ->■ a S pH{r,t) = ^ AiW(r - rp h), (103) 

i 

where as before, W is normalized, fW( r — r'; h)d 3 r' = 1, and having the property of finite support, 
W(r — r']h) —> 0, for |r — r'| > h, and the weight Ai should be chosen appropriately to minimize 
the difference between a(r,t) and aspn(r, t) everywhere. The above expression means that we are 
representing the continuous density as sum of finite number of unit distributions (kernel) carrying 
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the quantity A*. These unit density distributions are centered at the position r, (t). Due to the 
normalization of the kernel W, we have 


N 

/ a SPH {r : t)d 3 r = ^Ai, 


which should be the total value of the extensive quantity A of the system. 

For the application in hydrodynamics, we can use these parameters as the variational variables in 


Eq. (96) by representing the conserved density, for example the entropy density s, in the SPH form as 

N 

ssHp(r,t ) = y^ VjWjr ~ h ), (104) 


where the weight u^s are constant in time. Another extensive quantity, say A , is calculated as in Eq. 


(JT03j) with weights, A, = ( a/s) i Pi . In this scheme, the conservation of ssph is automatically satisfied 

j sph{?A) = y ^ijViW(Y ~r j(t)). (105) 


by defining 


The effective Lagrangian Eq. (96), written in the SPH representation, is 

' E' 


L S ph ({r*, rj) = - ^ v i (e/s) i = - 


7 


(106) 


where E t = (e*/s*) i is the “rest energy” of the i-th “particle” which carries the amount of entropy z/j. 
The equations of motion are obtained from the variational procedure as 


d f Pi + 

— I Vi -7i r, 


dt 
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ViVj 


P l 


Pi 


—— ——— 
a *2 ^ *2 
L 7 *7 . 


VjVF(r j - vf h) = 0 . 


(107) 


One of the advantages of the SPH approach is that we do not introduce numerical spatial derivatives 
using a difference method. In fact, in SPH spatial derivatives can be calculated analytically in terms of 
the kernel W. Another nice feature of SPH is its flexibility since it can be used in problems with initial 
conditions without any symmetry. Furthermore, due to the Lagrangian nature of the method, it is 
suitable for explosive processes such as those found in relativistic heavy ion collisions and it is also very 
easy to construct the freeze-out surface in this method. In addition, the variational approach guarantees 


that the SPH equations (107) give the optimal equation of motions for the parameters (r^t)} within 
the given total number of “particles”. In this approach, no numerical instabilities hardly occur since 
the whole system is a Lagrangian system and the total energy of the system is in principle conserved. 


Chiral Field and Hydrodynamics The variational formulation [ 383] can be used to obtain the 
equations of motion of a relativistic plasma of quarks interacting with the mass generating chiral mean 
held at finite temperature |384H392] . The effective Lagrangian density for the four-component isovector 
chiral held, (ft = (cr, 7f), where a is a scalar held and tt 1 are pseudoscalar helds playing the role of the 
pions, in the presence of the quark thermal bath which is at rest is expressed as 


4// = 


(108) 


where is the thermodynamic potential. On the other hand, the effective Lagrangian in Eulerian form 
with constraints from the conserved quantities is 


C 


(fluid) 

eff 


= -e(n, s ) - nu^d^X - su^d^C + ^w (u^u^ - 1) , 


(109) 
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where n and s are the conserved proper number and entropy densities and A, Q and w are Lagrangian 
multipliers. It can be shown that the Lagrangian density evaluated in the proper comoving frame of 
the fluid is equal to [383] 

jrUJ™d) _ _ e ( nj s) + im + Ts — p = —Q. (HO) 

Our coupled system is then described by an effective Lagrangian 

£(4+fluid) = _ e ( Wj S} (j)) - nu^d^X - su^d^C + ^ w (tt'X - 

i). (in) 

and the corresponding equation of motions are 


□ ( I> = —R, d u T^ uid) = Rd^ , 

(112) 

where R is given as 

R = ' + Pi 4>) , 

(113) 

with 

f d?k 1 /EM , , , . , 

9 Uq J (27r) 3 + \ + 9q) 

(114) 


is a scalar density for quarks. Here u q stands for the color-spin-isospin degeneracy of the quarks, 
E k (4>) = ( k 2 + and m q {4>) = (g 2 0 2 ) 1,/2 = + v ? 2 ) 1 / 2 plays the role of an effective mass for 

the quarks. 


Viscosity and Variational Approach When the fluid is not ideal, the above mentioned varia¬ 
tional method requires the introduction of some other components to take into account the violation 
of time-reversal symmetry since the standard variational method is always time-reversible. One way to 
accomplish this is to introduce a Rayleigh dissipation function and extend the variation of the internal 
energy to include dissipative effects. In the context of general relativity, the Lagrange approach for 
viscous hydrodynamics has been developed [57511353] . There, heat conduction and multiple component 
fluids are discussed. The SPH approach can also be extended to include viscosities by introducing the 
time variation of the weight factor {zy} and equations of motion for viscous tensors [ 162 . 394 , 1395 ]. 

Another way to introduce dissipative phenomena in the variational formalism was developed by K. 
Yasue by replacing classical variables by stochastic ones [ 355] in the so-called Stochastic Variational 
Method (SVM). It has been shown that the non-relativistic NSF equation can be derived from the 
SVM 13971 and this method also opens up an interesting possibility for the interpretation of the origin 
of quantum mechanics [397], However, unfortunately its relativistic form is not yet known and, in 
the pragmatic side, an application of the SVM in solving viscous hydrodynamics equations requires 
technically non-trivial procedures. 

6.3 Known Analytic Solutions 

Several analytic solutions of relativistic hydrodynamics are known. Naturally, they describe some spe¬ 
cific situations but they often offer important physics insights and in some cases they serve as a powerful 
tool to study realistic situations even at the quantitative level. Furthermore, analytic solutions can be 
used to check numerical codes, which is of fundamental importance in numerical calculations of rela¬ 
tivistic hydrodynamics since they employ many subtle technical details to avoid numerical instabilities. 

The first analytic solution of ideal relativistic hydrodynamics in (1+1)D was given by Khalatnikov 
[23] for the massless ideal gas initially at rest and confined within a finite spatial interval with thickness, 
say 2£. This initial condition was used in Landau’s model for the multiple pion production [22, 23] , 
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By changing the independent variables (t,z) to the hydrodynamic variables (a, (), where a is the 
fluid rapidity, and ( = In (T/T 0 ), the light-cone variables, (t ± z) ft are expressed in terms of integrals 
of modified Bessel function Iq containing dimensionless variables a and (. Therefore, to obtain the 
distributions of thermodynamic quantities and the fluid velocity as function of space and time, we 
still need to invert the resulting coupled transcendental equations. In spite of this complexity, the 
Khalatnikov’s analytical solution is useful to extract physical insights of the flow. Landau predicted 
that the rapidity distribution of the final produced particles was approximately Gaussian with a width 
that was a function of the incident energy. Further studies and applications of this model have been 
extensively done, including the generalization of the EoS from a gas of massless particles to an arbitrary, 
but constant speed of sound (? s [3981 - 1402] . The observed Gaussian rapidity distribution of charged 
particles in A-A collisions has been discussed within the scope of the Landau model [40211405] . For a 
more recent review on the Landau model and the Khalatnikov solution see [406] . As an application of 
the Landau model in the case of constituent quarks, the system size dependence of the multiplicity and 
the rapidity distribution are discussed in Ref. [407] . 

The Landau model and its analytic solution become very simple in the high energy asymptotic 
regime and at very high energy the central rapidity regime becomes boost invariant [408H410] . In this 
case, the rapidity of the fluid element is given exactly by 


1 t + z 

y — ri = - In-. 

y 1 2 t-z 


(115) 


As mentioned in Sec. 4.1 Bjorken argued physically that such a situation should occur in relativistic 


heavy ion collisions at ultra-relativistic energies and derived a simple formula to estimate the energy 
density in the initial state [116] . Although now we know that the real situation does not exactly satisfy 
those scaling conditions, this has been used in many (2+l)D calculations in the central rapidity domain. 
It is also emphasized in [ 411] that, within realistic values of parameters in A-A collisions, the time scale 
for which the longitudinal flow approaches the Bjorken scaling limit starting from the full-stopping 
initial condition is too large even at LHC energies, which suggests the necessity of some longitudinal 
flow profile already in the initial state. 

In order to take into account the non-boost invariant nature of heavy ion collisions, Bialas, Janik, 


and Peschanski m generalized Eq. ( |115[ ) as 2 y = In f + (t + z) — In /_ (t — z ), where f± are functions 
to be determined by the hydrodynamic equations. These functions obey an identical simple first order 
differential equation, and from this they have obtained a one-parameter family of exact solutions for 
relativistic hydrodynamics in ( 1 + 1 ) dimensions, which interpolates between the two limiting solutions 
of a Hwa-Bjorken boost invariant case and the Landau-Khalatnikov full-stopping initial condition [412]. 

The Budapest group (Csorgo, Nagy, Csanad) found a simple family of accelerating solutions in ideal 
relativistic hydrodynamics that are exact, explicit and analytic [41311416] . by assuming that the solution 
depends only on the Rindler coordinates, t = r cosh 77 , r = rsinh / 7 . Their solutions are expressed as 


v = tanh(A? 7 ), 


P_ 

Po 


?) 


\d(K-\-l) / K, 


cosh. Md -‘> (§) 


(116) 


where v (r, 77 ) and p (r, 77 ) are the velocity field and pressure, respectively, and d is the dimension of the 
space. It is important to note that the coordinate r for d > 1 is the radius of a d -dimensional sphere. 
The constants parameters A, k and (j)\ specify the type of physical solutions (in particular is the 
velocity of sound squared, specifying the EoS as before), but not all of these constants are independent 
and they are constrained for given value of d. Depending on the constraint, some of these solutions 
reduce to already known solutions [ 1161 14081 141711420] while others are new and accelerating ones. 
The accelerating solution can be used, for example, to estimate the initial energy density from the 
experimental observed rapidity distribution, and they showed that the Bjorken formula significantly 
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underestimates the real value at RHIC A-A collisions since it neglects the effects of work done by 
acceleration. For a more detailed discussion we refer the reader to Refs. [ 4 15 L I416J . 

Pratt developed a method to solve the equation of motion of relativistic hydrodynamics in the co¬ 
moving (Lagrangian) coordinates and obtained an explicit solution for the (1+1)D dimensional case 
with an arbitrary value of the speed of sound (the parameter n above) for the initial Gaussian density 
profile. He analyzed the longitudinal acceleration of this model, in particular with respect to the HBT 
measurements of longitudinal size [ 4211] . 

All of above mentioned analytic solutions are essentially (1+1)D solutions due to symmetry even if 
the space has d spatial dimensions. The first (3+l)D simple, explicit analytic solution for cylindrically 
symmetric case was found by Biro m, where the transverse expansion of the fluid involved the 
equation of motion corresponding to a mixed phase with a first order phase transition, P = const. The 
solution was also extended to include non-zero initial velocity profiles mg. In this case the solution is 
still analytic but not anymore explicit. Sinyukov and Karpenko [ 422 ] studied the same EoS to extend 
the ellipsoidal case to study elliptic flow in ideal hydrodynamics. 

Khalatnikov’s method to obtain transverse flow with longitudinal boost invariance was carried out 
by Peschanski and Saridakis [423j assuming that the transverse flow is quasi stationary and driven by the 
cooling of the temperature due to the longitudinal expansion. Under this condition, the transverse flow 
pattern was calculated analytically by a similar method used in the Khalatnikov solution. Identifying 
the entropy flow as particle flow, they studied the relation between the initial eccentricity and the final 
elliptic flow coefficient. 

A different class of cylindrically symmetric analytic solutions was constructed by Lin and Liao [ 424 ] 
by embedding a transversal Hubble flow in the (3+l)D equation with an EoS of the type e = up and 
then constructing the longitudinal flow equation. In addition to the 2D and 3D pure Hubble solutions, 
they obtained a set of non-trivial explicit solutions. 

The method based on self-similar flow initially introduced to find solutions for non-relativistic cases 
known as Buda-Lund approach [ 1701 I425H434] has shown to be very powerful in studying the flow 
dynamics of non-relativistic hydrodynamics. Recently this approach was extended to include (rigid) 
rotation, and explicit solutions with rotation were found [ 435 , I436j . In this line, a class of analytic 
solutions for truly 3-dimensional solutions for ideal hydrodynamics were obtained by generalizing the 
Bunda-Lund self-similar flow method to the relativistic cases [ 419 . 420 ]. 

The key point of this approach is that the continuity equation in d-dimensions can be solved for 
some specific flow profile, p (r, t) = g (t) f (r 7 A _2 r), where A is a time dependent diagonal d x d matrix 
A = diag (Ad, X d ) , and / (s) is an arbitrary smooth function of the scaling variable s = r T A~ 2 r with 
normalization J 0 °° s d ~ 1 ds f (s) = 1. The velocity profile is then given by v=A~ 1 A (t) r. It is easy to 
verify that these ansatze for p and v satisfy the continuity equation for any {Ad (t) ,.., X d (t)}. Taking 
the EoS as that of an ideal gas of massive quanta, 

£ — mn + £ in , p = nT , (117) 

accounts for the presence of mass flow. The solutions given in Refs. [ 419 1 1420] have a simple, explicit 
analytic form and correspond to truly (3+l)D cases that contain the Hwa-Bjorken solution as a special 
case. In this scenario, it was shown that the final observables are directly related to the initial distribu¬ 
tion function / and solutions of the Buda-Lund type predict a universal feature of the observables. Such 
aspects were discussed in Ref. [437 ] . Recently, the Budapest group further extended these solutions to 
incorporate general, realistic equations of state [ 438 ] and also multipolarity of the initial distribution 
to analyze the experimental data of higher order flow coefficients [439] . 

These studies show that the choice of certain symmetry aspects of the solution play a fundamental 
role in constructing analytic solutions. In the limit of conformal symmetry, a class of analytic solutions 
of the transverse flow with longitudinally boost invariant background was first shown by Ref. j 1 IQj . 
The conformal symmetry assumption also allowed to include first order viscous corrections [ 440] . This 
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approach has been developed offering a more general theoretical scheme to study a variety of physical 
situations, including angular momentum and anisotropic flow in a perturbative form [431] . 

The generalization of Gubser’s work to Israel-Stewart theory (IS) was carried out in Marrochio et al. 
[4421) and the first analytic (and semi-analytic) viscous solutions were obtained for the boost invariant 
transversely expanding flow. These IS solutions [442] have well defined temperature and shear stress 
tensor profiles and do not display negative temperatures as in the case of the first order theory originally 
considered in Ref. [441] . In fact, in Ref. [ 442] a comparison between the first order case and IS was 
performed and the importance of the relaxation time was discussed. More recently, Hatta, Noronha, and 
Xiao and collaborators in their series of works [44311446] explored various aspects of conformal symmetry, 
finding many new analytic solutions for viscous relativistic hydrodynamics in different situations and 
their properties were discussed. In particular, in Ref. [444] . they showed that several new solutions for 
second order viscous hydrodynamics can be obtained using different background geometries, which also 
include vorticity effects. Furthermore, Ref. [ 445 ) calculated analytically the harmonic flow parameters 
as perturbations around radial flow and obtained explicitly the n dependence of v n and discussed also 
the effects of shear viscosity (a semi-phenomenological discussion on bulk viscosity was also included). 

As we have seen, these solutions have played an important role in the sense that they offer many 
checking points of the physical assumptions assumed in hydrodynamic modeling. In analytic approaches, 
obviously the real physical event by event fluctuating observables are not incorporated by definition. 
However, some well-chosen observables reflect the symmetry of the system so that the analytic solutions 
which meet these symmetry requirements can be discussed in the framework of analytic approaches. 
For such purposes, explicit and analytic solutions from which physical observables can be calculated 
also analytically are extremely useful. For example, the Buda-Lund type of solutions clearly exhibit the 
observed scaling laws [438] . such as the transverse energy dependence of HBT radii, mass dependences 
of the effective temperature of single particle spectra, and etc. They also serve to study the fundamental 
problems intrinsic to the hydrodynamic approach, such as the onset of instabilities in relativistic viscous 
hydrodynamics for a given equation of state and transport coefficients. 


6.4 Miscellaneous Aspects of Relativistic Hydro Modeling 


Hydrodynamics represents the macroscopic motion of a fluid close to local thermal equilibrium. If the 
system is far from equilibrium, it is not easy to describe the dynamics correctly only with macroscopic 
dynamical variables. In fact, in numerical simulations with large values of 77 /s or (/s the corresponding 
viscous correction to the single particle distribution function 5f becomes easily as large as the thermal 
equilibrium function f eq . This means that the basic assumption behind the kinetic derivations of 
dissipative hydrodynamics does not hold in these extreme cases. However, in some specific situations we 
may expect that the macroscopic degrees of freedom can still cover the principal features of the dynamics 
of the system without resorting to a full microscopic description of the non-equilibrium dynamics. 

This is the case of the very early stage of the hydrodynamical evolution where the there is still a clear 


distinction between the transverse and longitudinal directions as we mentioned before in Sec. 4.1 and 
also will be discussed in Sec. [ 8 ] In such a situation, the single particle distribution in momentum space 
will not be isotropic in contrast to the case of LTE. The extreme case where the longitudinal direction 
is free-streaming while the transverse plane is thermally equilibrated was first proposed in Ref. [447J. 
See also Ref. [448] . Later, a theoretical scheme of hydrodynamics which takes into account such an 
anisotropy in the momentum space including dissipative effects has been developed in Refs. [449H45T] . 
For comprehensive reviews see Refs. [452 , 453 ] and the recent review of relativistic hydrodynamics by 
Jeon and Heinz where detailed discussions on this approach [454] are given. Longitudinal/transverse 
difference in momentum space can also occur when the incident flow of the colliding matter is incor¬ 
porated and such a situation is generally related to the discussion of rapidity dynamics of baryon rich 
matter. To handle such a situation, multicomponent fluid models have been studied [4551J457] . These 
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models are considered to be appropriate in the description of the rapidity dependence of the direct flow 
parameter V\ in the energy regime relevant for NICA and FAIR experiments. 

A quite different view of the origin of the deviation of the single particle spectrum from that of 
Boltzmann comes from the so-called non-extensive statistics. It has been long known that the single 
particle spectra of particles in high energy collisions decay exponentially as function of the transverse 
energy for lower values of y 7 s but, at higher values of they have a power law tail. In view of QCD, 
such power law tail of hard components can be well understood in terms of perturbative processes, and 
the soft exponential part can be understood in terms of canonical statistic effects. Phenomenologically, 
both behaviors can be interpolated by the so-called Levy-Tsallis distribution [4581 - 1463] . Tsallis derived 
this distribution postulating a non-extensive nature of the entropy [464]. Some authors have discussed 
the possibility to incorporate such a formalism in the construction of dissipative hydrodynamics [ 4651 
1466] . Osada argued that such a situation can be considered as a consequence of the existence of long- 
range correlations which introduce non-extensive behavior in the local thermodynamic relation [41)7] 
and also proposed a generalized matching condition including such a situation with the presence of 
particle and anti-particle distribution functions [468 j. 

An important point that has been less extensively discussed is the effect of angular momentum 
and the possible emergence of instabilities in fluids with low viscosity. In the usual hydrodynamic 
description of relativistic heavy ion collisions, this subject has attracted little attention since it requires 
a high resolution calculation with a very small numerical viscosity. However, if we look for genuine 
hydrodynamic signals, the nonlinear response which leads to instabilities (for example, Kelvin-Helmholtz 
instability) may offer important additional sets of observables as advocated in Ref. [469] . The study of 
relativistic hydrodynamics with spin was developed by Becattini [ 470] . The use of the effects of rotation 
and associated phenomena such as vortices and turbulences as new tools for precise understanding 
of hydrodynamic modes has recently been discussed in Ref. m- Recently, exact solutions of the 
relativistic Boltzmann equation in the relaxation time approximation (RTA) form have been constructed 
and used to derive dissipative hydrodynamic equations, including anisotropic hydrodynamics [472] , 
Later, solutions of the RTA Boltzmann equation corresponding to Gubser flow were constructed and 
the dissipative hydrodynamics was derived [47311474] , These interesting subjects are quickly developing 
and are unfortunately beyond the scope of the present review (see, for example, Refs. [32511475] ). 


7 Dynamical Effects of Coarse-Graining in Hydrodynamics 


So far we have seen that, in spite of the great success of the hydrodynamic description of relativistic 
heavy ion collisions and the general picture of the Little Bang being established as some kind of “stan¬ 
dard model” of heavy ion collisions, the proper physical foundation of the validity of hydrodynamics 
in relativistic heavy ion collisions is yet far from being established. The main quantitative question 
involves the validity of the LTE on an EbyE basis and the size of the coarse-graining (C-G) scale as¬ 
sociated with it. As we have mentioned, the physical observables at the moment are yet not sufficient 
to determine uniquely the hydrodynamic picture. Furthermore, as suggested in Sec. 6 .2[ hydrodynam¬ 
ics may also correspond to an ensemble averaged dynamics, in accordance with the observables from 
the experimental measurements, so it is crucial to understand how precisely one is able to probe the 
system’s profile. If the observables investigated do not have enough resolution of the macroscopic dy¬ 
namics, then the properties of the matter deduced from them can be just an effective description, e.g. 
corresponding to an effective EoS and transport coefficients. Therefore, one needs to be careful with the 
assumptions required for the validity of the hydrodynamic description and whether they are supported 
by experimental observables. A precise analysis of the profile of the system’s dynamics may provide 
essential information about the properties of the matter as well as the initial condition of the collision 
(cf. discussions in Ref. [476] ). However, using only the available data, it is a difficult task to reconstruct 
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the event profile evolution. 

On the other hand, there exist several transport models which exhibit collective behavior simi¬ 
lar to the matter created in heavy ion collisions. In particular, the PHSD model developed by the 
FIAS/Frankfurt and Giesen Groups can reproduce well the observed collective signals (see for example, 
Refs. [1401114211477H479] ). but the LTE is not obviously being explicitly implemented in this approach. 
Therefore, it is interesting to study the hydrodynamic aspects of such a model by introducing the hy¬ 
drodynamic variables through the C-G procedure described in the Sec. |6.2[ In this way, we expect to 


test the correlation between the observed collective flow parameters and the concept of LTE used in 
the hydrodynamic approach. 

In the PHSD model, the transport dynamics are described in terms of quasi-particle modes of the 
Kadanoff-Baym equation [480] . and the effects of many-body interactions are taken into account in 
the form of a mean held, which in turn, generates the mass of the quasi-particles. Therefore, this 
model is particularly suitable for our purposes since it represents the quantum mechanical microscopic 
transport dynamics in terms of effective (resummed) strongly interacting dynamical quasi-particles and 
off-shell hadrons as degrees of freedom [14011143] . In this section, we report shortly on the hydrodynamic 
representation of the PHSD dynamics and see how much the LTE is attained dynamically. 


7.1 Implementation of Coarse-Graining in PHSD 

A single event in the PHSD simulation provides the evolution of coordinates and momenta corresponding 
to the quasi-particles. Therefore, it is straightforward to apply the C-G procedure similar to Sec. |6.2 
to obtain the macroscopic energy and momentum tensor, 


T w (x) = 


P\ 


,o 


r i(t);h ), 


where IE is a normalized smoothing function (here we take a Gaussian with width h) and the sum is 


taken over all the quasi-particles of the event. However, as mentioned in Sec. 6.2, we can also include 
the different “events” in this sum. 

Roughly speaking, an event of the PHSD simulation is a kind of particle observation for a given 
quantum wave function of the system. However, a single event does not carry information about 
the interaction among them. Therefore, we need to introduce an ensemble of such events (referred as 
“parallel events”) to take into account the dynamics of the wave-function which contains the information 
of many-body interactions. Such information is determined by the particles of a number of parallel 
events starting from the same initial condition. This number, referred here as NUM, represents the 
number of ensembles of parallel events to generate the mean-held effects. 

By increasing NUM, the effect of the potentials shows up in the dynamics of the system’s evolution. 
On the other hand, if we take numbers that are too large, the corresponding energy density profile looses 
its particle nature, that is, we loose the resolution in the classical sense. Thus, NUM also constitutes a 
parameter of the C-G procedure in the present approach. 

We have performed several tests using a simulated sample of events produced with PHSD in order 
to check the hydrodynamic properties extracted from the C-G procedure described above. Here, we 


summarize the results presented in Ref. [481] . The plots in Fig. 12 show the longitudinal profiles of the 
components of the particle velocity for two snapshots of the evolution (the open circles represent the 
mean values of the distributions) for one single PHSD event of a Au-Au collision at y/sjsrN = 200 GeV. 

We tested two types of kernel functions: (a) a simple rectangular box in which W=l/(AxAyAz) if 
the particle is inside the box and IU=0 if it is outside, with Ax=Ay being the transverse lengths and 
A z{t) the longitudinal length of the box; and (b) a 3D gaussian function with characteristic lengths 
defined by the widths h t in the transverse direction and h z {t) in the longitudinal direction. The left plot 
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Figure 12: (color online) Snapshots of the longitudinal profile of the particle velocity in PHSD for two 
different time-steps of the evolution for a central Au-Au collision at yj snn = 200 GeV. The color scale 
in each panel represents the number of particles. 


in Fig. 13 shows the dependence of the local energy density as a function of time for the rectangular 
box (red curve) and the Gaussian kernel (black curve). The energy-momentum tensor in this case 
was calculated in the neighborhood of the origin (x=0, y= 0 , z—rj—0 ), considering a rapidity window 
A?/=0.1. We tried to choose the sizes of each kernel to be roughly of the same order though they are not 
exactly equivalent. As expected, both cases presented similar behavior with the Gaussian kernel being 
smoother and more strongly correlated. The plot on the right hand side of Fig. 13 shows the dependence 
of the ratios of the longitudinal (Pl) and transverse (Pt) components of the pressure with respect to the 
proper energy density e as a function of the time evolution of the system. In this case we computed the 




Figure 13: (color online) Dependence of the local energy density e with the time evolution for two 
different types of kernel functions (left plot); dependence of the longitudinal and transverse pressure 
components with the system evolution considering different characteristic lengths in the longitudinal 
direction (right plot). 

T ,w in the neighborhood of the point (x, y, z)—(l, 0, 0) fm (avoiding the divergence at the origin) and 
four different sizes of the longitudinal length of the Gaussian kernel were tested, which are displayed in 
the legend of the plot by Ar] and corresponding to the range of — cr to +cr around the Gaussian peak. 
The transverse length was fixed as 0.3 fm (which corresponds to la of the Gaussian kernel). For all 
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configurations tested we noticed that the transverse components of the pressure always start very low 
and then increase as the system evolves, while the longitudinal component always starts very high and 
then decrease. For A?/=0.3 (dotted lines) the pressure components become isotropic and last for an 
interval of about 2 — 3 fm/c, starting at f~1.0 frn/c. However, with a further increase of the longitudinal 
length of the Gaussian kernel, the isotropy of the pressure components is not maintained. Actually, 
the large size of A// introduces an artifact due to the mixing the collective longitudinal expansion. On 
the other hand, in the infinitesimal limit of A r], EbyE fluctuations become too large. Different sizes for 
the transverse length of the kernel function were also tested without producing any important changes 
in the observed behavior. Finally, we varied the parameter NUM (1, 2, 5, 10, 20 and 30) in order to 
check the effects due to the mean-field potential and we found that the general behavior converges for 
NUM>5, with the fluctuations being very large otherwise. 


7.2 Collectivity in the evolution of PHSD events 

In order to check whether the initial spatial anisotropy in the PHSD events is converted into a final mo¬ 
mentum anisotropy as one would expect in a scenario with collectivity, we evaluated the time evolution 
of the eccentricities e n of the particle distribution in the transverse plane as well as the flow coefficients 
v n computed with respect to the corresponding eccentricity phase angle: 


(r n cos (n [0 — <F n ])) 
£n = 0F 

v n = (cos (n [-0 - 'F n ])), 


<F n = — arctan 
n 

d'n = f hn + 7 r/n, 


(r n sin (n0)) 
(r n cos (n0)) 


(118) 

(119) 


where r = \Jx 2 + y 2 is the distance of a given particle to the center of mass of the distribution in the 
transverse plane, 0 = arctan (y/x) is the angle of its position vector with respect to the ir-axis, and 
0 = arctan ( p y /p x ) is the angle of its momentum vector with respect to the x-axis. 

The left plot in Fig. 14 shows the behavior of the coefficients £ 2 and V 2 as a function of the evolution 
of the system for a single PHSD event. The plot on the right shows the distribution of the final V 2 




Figure 14: Time evolution of the spatial eccentricity e 2 and the elliptic flow coefficient v 2 for a single 
event (left plot); the distribution of e 2 vs v 2 (right plot) for many events. The red line represents the 
average (u 2 ). 

versus the initial e 2 coefficients for a sample of many events. The positive correlation observed is an 
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indication of collective evolution of the system. Moreover, the relevant transfer of the anisotropy from 
the initial to the final states occurs during the initial (partonic) stages. 

Higher harmonics show similar behavior for central and intermediate cases, but the “quality” of the 
transfer becomes poor for the most peripheral cases. Very similar results have been observed in ideal 
hydrodynamic calculations as discussed in Refs. [21711482] . which demonstrates that even if complete 
LTE is not present, the usual observables that signal the system’s collective behavior are similar and 
not sensitive to deviations from equilibrium. Therefore, the development of new observables that are 
sensitive to deviations from equilibrium of the matter during the time evolution is of great interest in 
order to understand the collective properties of heavy ion collisions and the evolution of collectivity. 


8 Discussion and Future Problems 


The hydrodynamic picture has shown to be very successful in describing many aspects of the dynamics 
of the matter produced in relativistic heavy ion collisions. Such a picture is consistent with the Little 
Bang picture of the Quark-Gluon Plasma. In this work we have reviewed, from a pedagogical point 
of view, the basic structure and physical meaning of the relativistic hydrodynamic approach with a 
short historical overview. We mention also the most recent developments based on the LHC and RHIC 
observations. In particular, we emphasized the meaning of coarse-graining implicitly contained in the 
hydrodynamic description. Although the hydrodynamic description achieved an overwhelming success 
in reproducing the experimental data, there still remain some important questions associated with the 
potential existence of an unknown C-G scale, which should be clarified. 

For example, C-G is often loosely quantified in terms of the Knudsen number in kinetic derivations 
but reality is much more subtle because the kinetic approach itself contains a kind of C-G from the 
beginning. At first sight, it seems that the Boltzmann approach is somewhat conflicting with the idea of 
a strongly interacting and dense sQGP fluid because of the assumption of dilute gas, which is the starting 
point of the Stossahlansatz. In spite of this, the Boltzmann equation is considered to be very useful 
to discuss relativistic hydrodynamics in particular to study the structure of relativistically covariant 
dissipative terms. This is so because, due to the assumptions of local binary collisions and the absence of 
memory effects, relativistic covariance is easy to maintain. In the presence of non-local interactions and 
correlations due to three (or more) body collisions, it is quite difficult to satisfy relativistic covariance, 
as is known in the traditional cascade calculations [55]. In any case, in the derivation of hydrodynamics 
from a microscopic theory, the extraction of macroscopic variables and the associated C-G of space-time 
scales are the fundamental procedures. In addition, to close the system of equations in these macroscopic 
variables, we introduce the thermodynamic/statistical properties (usually) implicitly depending on C- 
G. Unfortunately, there is no established method to perform a C-G procedure which is applicable to 
any general situation in relativistic heavy ion collisions, and even a general method may not exist at 
all. In fact, the C-G procedure is not known even for non-relativistic systems far from equilibrium. 


Several approaches discussed in Sec. |6.1| are examples of these C-G procedures in different situations. 
In this sense, the formulation of relativistic hydrodynamics within the context of applications in the 
description of relativistic heavy ion collisions is not yet completely established and further investigations 
are necessary. 


As discussed in Sec. (L2 due to C-G one hydrodynamic event represents an extremely large number of 
microscopically distinct physical events. The larger the C-G scale is, the easier a macroscopic description 
such as hydrodynamics becomes applicable. However if a much larger C-G scale is used, we loose the 
required resolution in the space-time recognition. In fact, we cannot observe inhomogeneities with 
smaller wavelength than the C-G scale. This affects directly the class of observables that the model can 
describe. 

Even though, we do not know a priori what kind of observables are insensitive to the C-G scale. As 
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an example, we recall that the success of hydrodynamics was first established by the analyzes of the 
elliptic flow coefficient using smooth initial conditions (s-IC) and, subsequently, EbyE analyzes with 
fluctuating initial conditions (f-IC) confirmed this hydrodynamic interpretation for v 2 . But this seems to 
be somewhat amazing because the hydrodynamic equation couples the fluid flow and the corresponding 
driving force in a very complicated nonlinear form. If elliptic flow is due to genuine non-linear signature 
aspects of hydrodynamics, one might expect that the two analyses should give appreciable differences. 
For example, if we consider an observable O as functional of the initial condition IC, in general we 
expect for a non-linear response, 

(o [IC\) + O [</C)l, 


due to the hydrodynamic evolution, where () represents the average on different IC’s. In the equation 
above, on the right-hand side we identify s-IC as an average over all f-IC fixing a parameter such as 
centrality. 

I 11 fact, some differences between the left and the right-hand sides have been observed for some 
observables such as the p t dependence of u 2 , particle spectra, and HBT radii, but there is almost no 
difference in the integrated v 2 , except for very central collisions. These differences can easily be absorbed 
in the systematic uncertainties of the hydrodynamic modeling. This fact can be understood since for 
V 2 the effects from the initial pressure gradient, which is related directly to e 2 , is overwhelming over 
the small inhomogeneity in the initial conditions. This picture was quantitatively demonstrated in Ref. 
[ 483 ] using Monte-Carlo Glauber type initial conditions. As shown in Fig. 15 we first observe that these 
plots indicate that the average elliptic flow is very much insensitive to the fluctuations in the initial 
conditions since the width of the distribution of events around the average is very narrow for a given 
value of e 2 . Second, the event-averaged v 2 has a very keen linear dependence on e 2 - Therefore, any 
nonlinear aspect of genuine nonlinear hydrodynamic response is not contained in v 2 . In other words, 
the information of the detailed IC through v 2 is limited up to the resolution specified only by the initial 
eccentricity e 2l and any other inhomogeneity is hindered. In fact, if we use the integrated v 2 as a filter 
for initial condition selection, the corresponding C-G is too large and the hydrodynamic information 
obtained only from this observable is rather imprecise. The fact that v 2 is almost linear in e 2 and does 
not show significant EbyE fluctuations indicate strongly that v 2 is almost uniquely determined by the 
initial geometry. As is pointed out in Ref. [484]. several puzzling “scaling properties” are observed, in 


particular the quasi-independence of v 2 (pt )/^2 with energy and system size, and the very good scaling 
of (v 2 ) rv_/ (Pt) (1 /S)(dN/dy). It is interesting to investigate whether these observations can be 

understood in hydrodynamics, and the role of C-G 



Figure 15: Correlation v 2 vs e 2 in fluctuating initial conditions. Figure^ taken from Ref. [483]. 

The higher order harmonics are more sensitive to the effects of inhomogeneities, as seen from the 
fluctuation measure in scattered plots of v n vs e n [483]. Of course, a more detailed structure of v 2 may 

“Reprinted with permission from H. Niemi et al., Event-by-event distributions of azimuthal asymmetries in ultra- 
relativistic heavy-ion collisions, Phys. Rev. C87 (2013) 054901. Copyright © 2013 by the American Physical Society. 














also give a larger fluctuation, for example at large p?, especially if the initial state dynamics is more 
peculiar than that of the Monte-Carlo Glauber one generating very sharp hot spots. 

As mentioned in Sec|5j in LHC experiments the anisotropic flow parameters can be determined in 
EbyE basis and their distributions can be measured, which were not possible in RHIC experiments. Due 
to the strong correlations between {u n } and {e n } in hydrodynamic calculations, the LHC data on flow 
parameters provide unique insight into the initial state of the matter produced in relativistic heavy ion 
collisions. In fact, one can classify the event geometry class more accurately than only using centrality. 
In particular, not only the distributions, but also the correlations between v n and v m or their higher ones 
and event-plane correlations are considered to be useful to map the final state to the initial condition 
quite precisely [ 160 . [377 ;. 11831 i485j| . The opportunity of observing collective signals for one single event 
opens several remarkable possibilities associated with genuine nonlinear hydrodynamic behavior. In 
particular, one needs to construct observables which are sensitive to the traditional genuine hydro¬ 
signals, such as shock waves |486b488] and Raylcigh-Helmholtz instabilities [69] in terms of correlations 
of flow parameters in one single event. 

Most of the discussions about the hydrodynamic behavior of the matter produced at RHIC and LHC 
until now are concentrated on the mid-rapidity region. In this case, the boost-invariant aspects have been 
emphasized. This is consistent with the CGC-Glasma type initial conditions where the classical gluon 
field is first generated (boost invariant) and then is subsequently converted into partons. As mentioned 
in Sec. 4.1, the pre-equilibrium energy-momentum tensor possesses anisotropic pressure components 


and the transverse pressure is dominant over the longitudinal one. On the other hand, as shown in Fig. 
13 in the initial condition based on the PHSD model the longitudinal pressure in the corresponding 
stage is much larger than the transverse one. This difference comes from the fact that in the PHSD 
approach the energy-momentum tensor is composed of particle-like objects (partons and strings). The 
huge longitudinal pressure is generated by nearly free-streaming particles going in opposite directions. 
This picture is not boost invariant at least in the very early stage where the usual hydrodynamic picture 
is not valid in the sense that the system is not in LTE. However, as discussed in Ref. [ 452] such a large 
value of the longitudinal pressure decreases very quickly due to the rapid longitudinal expansion and the 
transverse components dominate over the longitudinal one at f ~ 0.2 fm. At this time, the transverse 
components become more or less isotropic indicating that the anisotropic hydrodynamic picture may 
hold. On the other hand, in the CGC-Glasma picture, the energy-momentum tensor comes from the 
partons generated from a boost invariant classical gluon field. These partons have mainly transverse 
momenta. These apparently opposite behaviors represent a fundamental difference between the initial 
conditions for hydrodynamics constructed from transport and the CGC-Glasma pictures. It should be 
noted that the PHSD and CGC-Glasma approaches are based on QCD to obtain the dynamics of the 
quantum wave function of the collisions, which possesses the wave-particle duality. In our oppinion, it 
seems that the PHSD approach emphasizes the particle nature in this dynamics while the CGC-Glasma 
focuses on the wave nature (classical field). Although this point represents an interesting question, any 
clear observable which distinguishes qualitatively the difference between the two pictures has not yet 
been proposed to the best of the authors’ knowledge. 

In the above discussion, boost invariance plays a fundamental role to distinguish the two pictures. 
On the other hand, this refers only to the matter produced from the vacuum in the central rapidity 
region. In fact, thermal model analyzes show that the baryonic chemical potential is close to null at 
RHIC and LHC, which shows that the particles come dominantly form the excitation of the vacuum. 
However, for \y\ > 1, the rapidity distributions of charged particles are clearly deviated from the boost 
invariant form even at RHIC energies [ 489] . Furthermore, for the low energy heavy ion collisions such 
as BES at RHIC and the CBM experiment where one expects a large stopping power of colliding nuclei, 
the effects of finite baryonic chemical potential will not be negligible and the boost invariant postulate 
is not adequate. In this aspect, insight from analytic approaches such as those in Refs. [ 412.1415] will 
be useful as mentioned in Sec. 16.31 
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The determination of the flow parameters based on an EbyE analysis, even for photons and leptons 
in LHC experiments, also opens the possibility of investigating the time evolution of the collective 
profile realized in relativistic heavy ion collisions. As we have discussed, the proper formalism of 
relativistic hydrodynamics is still an open issue. One does not even know if a consistent framework 
exists that allows for the derivation of relativistic dissipative hydrodynamics for strongly interacting 
matter because, in this case, the concepts of the theory of relativity and thermodynamics require two 
different conflicting conditions: strict locality vs. thermodynamic limit. Studies of such observables 
related to the time evolution of the sytern will certainly offer a stringent check for the hydrodynamic 
picture and contribute to the determination of the transport properties of the fluid in question. 

In this review, we have explored several questions associated with relativistic hydrodynamics applied 
to relativistic heavy ion physics. In spite of the overwhelming successes of the hydrodynamic description 
of the experimental data from SPS-RHIC to LHC energies, we pointed out that there exist some 
disturbing questions from the conceptual point of view, mainly related to the role of coarse-graining 
and its intrinsic scale. There are some intriguing phenomena, for example, the observed flow patterns 
in p-p and p-A collisions and the associated ridge structure in two-particle correlations. If they really 
fit in the hydrodynamic picture, then the associated coarse-graining scale should really be small. If 
this is so, the realization of local thermal equilibrium occurs almost at the microscopic scale compared 
with the hadronic size, which would be rather amazing. In this sense, it is important to investigate in 
detail whether the thermodynamic relations and transport coefficients which quantitatively reproduce 
the data in such small systems remain exactly the same as those used in A-A collisions on an EbyE 
basis. 

Additionally, we pointed out that the boost-invariant approximation will break down when the 
nuclear stopping power becomes more effective. In this case, the incident baryonic flow becomes more 
dominant and we expect that some new features may appear compared to the dynamics of the matter 
directly created from the QCD vacuum. 

This field of physics is presently very quickly developing and new ideas and techniques are being 
actively proposed. We believe that the several questions raised in this paper will be crucial to improve 
our understanding of the physics of relativistic heavy ion collisions and QCD dynamics in extreme 
conditions. 
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